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ABSTRA.CT 



The development of a numerical procedure for the treat- 
ment of nonsimilar, unsteady, laminar boundary layers is 
presented. The solution is obtained from the laminar, iso- 
thermal, incompressible boundary layer equations employing 
a modification of the integral matrix procedure of Bartlett 
and Kendall. Solutions of example problems are presented 
for steady Blasius and Howarth flows, and for oscillating 
Blasius flow. Agreement with the known classical results is 
satisfactory and establishes the general feasibility of the 
method. Core storage requirements of 130,000 bytes allow 
consideration of as many as 25 nodal points across the 
boundary layer, 50 time increments per oscillation cycle 
and 50 streamwise stations. Solution of oscillating Blasius 
flow considering 8 nodal points and 16 time increments 
requires 13.49 seconds for one streamwise station utilizing 
IBM 360/67 time sharing capabilities. 
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TABLE OF SYMBOLS 





Rate of change of the mean free stream velocity 
with distance from the leading edge 
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Coefficients defined in the finite difference 
representation of streamwise derivatives 
(defined in equations (33) and (34) for two and 
three-point difference relations) 
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Stream function (defined by equation (14)) 
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Pressure 
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Local radius of the body in a meridian plane 
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or the leading edge 




Coefficients defined in the finite difference 
representation of time derivatives (defined in 
equations (37) and (38) for two and three-point 
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Time 
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Mean free stream velocity at the leading edge 
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u 


Velocity component parallel to the body surface 
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Velocity component normal to the body surface 




. Truncated series obtained in Taylor series 
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expansion of / f dn (defined by equation 
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Distance from the surface into the boundary 
layer measured normal to the surface 




. Truncated series obtained in Taylor series 
expansion of integrals involving nonsimilar 
terms (defined by equation (A7)) 
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a* 


Normalizing parameter (defined by equation 
(16) ) 
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Streamwise pressure gradient parameter (defined 
by equation (22)) 


k\-l 


Time difference between two time grid points 
denoted by the subscripts k and k-1 (defined 
by equation (39)) 




Logarithmic distance between two streamwise 
positions denoted by the subscripts i and j^-l 
(defined by equation (35)) 




Corrections for f.,f^,... during Newton-Raphson 
iterations 


6n 
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(defined by equation (30)) 


5t 


Time difference between two time grid points 
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Transformed streamwise coordinate (defined by 
equation (9)) 
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Density 


03 
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SUBSCRIPTS 
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Pertains to the boundary layer edge 


i 


Pertains to the i^^ nodal point in the boundary 
layer, starting with i=l at the surface 


k 


Pertains to the k^^ time grid in the oscillation 
cycle 


1 


Pertains to the streamwise position 
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m Pertains to the iteration during the 

Newton-Raphson iteration 

w Pertains to the wall or surface of the body 

SUPERSCRIPTS 

n Equal to unity for axisymmetric bodies and 

zero for two-dimensional bodies 

' (prime) Represents partial differentiation with respect 
to n 
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Coefficients of "linear" corrections in Newton- 
Raphson recurrence formulas corresponding to 
linear equations prior to any matrix reduction 
(defined by equation (61)) 

Coefficients of "linear" corrections in Newton- 
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nonlinear equations prior to any matrix reduc- 
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Newton-Raphson recurrence formulas correspond- 
ing to nonlinear equations prior to any matrix 
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Coefficients of "nonlinear" corrections in 
Newton-Raphson recurrence formulas correspond- 
ing to nonlinear equations after matrix reduc- 
tion (defined in equation (69)) 

Matrix of errors in Newton-Raphson recurrence 
formulas prior to any matrix reduction (defined 
in equation (60)) 

Matrix of errors in Newton-Raphson recurrence 
formulas corresponding to linear equations prior 
to any matrix reduction (defined by equation 
(61)) 
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[ela"^ 


Product of AlI ^ ^ -ElI (defined in equation 
(67)) - J ^ J 


i ENL j 


Matrix of errors in Newton-Raphson recurrence 
formulas corresponding to nonlinear equations 
prior to any matrix reduction (defined by 
equation (61)) 


Tenl ‘ 
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Matrix of errors in Newton-Raphson recurrence 
formulas corresponding to nonlinear equations 
after matrix reduction (defined by equation 
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[av] 


Matrix of corrections to primary variables in 
Newton-Raphson recurrence formulas (defined by 
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Matrix of corrections to "linear" primary 
variables in Newton-Raphson recurrence formulas 
(defined by equation (61)) 
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variables in Newton-Raphson recurrence formulas 
(defined by equation (61)) 
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I. 



INTRODUCTION 



A. GENERAL 

All aerodynamic flows are to some degree influenced by 
unsteadiness. Rotating stall in turbomachinery, aerodynamic 
stall flutter and hydrofoil flow fields are but a few 
practical problems in which the nonsteady fluid flow phe- 
nomena are of prime importance. Additionally, the recent 
emphasis on the development of rotary wing aircraft has 
created a new demand for knowledge of the behavior of 
boundary layers in unsteady flow. 

B. ANALYTICAL HISTORY 

Early analytical treatments by Stokes, Rayleigh and 
Schlichting (Ref. 1) considered only a special case of the 
general problem: that of unsteady viscous fluid flow with 
no mean flow or pressure gradient. Lighthill (Ref. 2) made 
the first significant investigation of small sinusoidal 
oscillations superimposed on a mean flow. Lighthill 's 
analysis was based on a small perturbation treatment of the 
boundary layer equations for sinusoidal flow, with only 
the first order terms retained. The determination of the 
existance of the "quasi-steady" regime at low frequencies, 
and the "shear wave" regime at higher frequencies was the 
principal result of this analysis. Glauert (Ref. 3) extended 
Lighthill 's work by considering the boundary layer in the 
vicinity of a stagnation point. Carrier and DiPrima (Ref. 4), 
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using a system of equations derived using a modified Oseen 
linearization, confirmed the existance of the "shear wave" 
solution and the phase advance of the wall shear determined 
by Lighthill. Nickerson (Ref. 5) initially retained per- 
turbation terms beyond the first order, thus widening the 
scope of Lighthill 's work. 

In addition to his analytical development, Nickerson 
made hot-wire measurements of the laminar boundary layer on 
a flat plate oscillating in a Blasius mean flow. His work 
partially confirmed both his own and Lighthill 's analyses, 
but mechanical difficulties precluded the completion of the 
experimental program. Hill and Stenning (Ref. 6) imposed 
sinusoidal oscillations directly on mean Blasius and Howarth 
flows and were able to confirm the analytical results of 
both Lighthill and Nickerson. Recently an integral analysis 
of the unsteady boundary layer was introduced by Koob and 
Abbott (Ref. 7) , who succeeded in demonstrating the method 
for the case of an accelerating flat plate. 

C. APPLICATION OF NUMERICAL METHODS 

With the advent of the high-speed digital computer, it 
has become possible to apply numerical methods to the 
analytical studies of viscous fluid flows. Fromm (Ref. 8) 
and Chorin (Ref. 9) introduced modified finite difference 
methods for the treatment of viscous flows which could 
feasibly be applied to the boundary layer problem. Smith 
and Clutter (Ref. 10) and Kleinstein and Nabi (Refs. 11 and 
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12) applied finite difference methods to the steady boundary 
layer equations. Successive approximation, linearization 
and integral methods have also been applied to the boundary 
layer equations (Refs. 13, 14 and 15). A major drawback of 
most of the methods mentioned is that in order to reach a 
solution to a given problem large computing times are nec- 
essary. 

The introduction of the Boundary Layer Integral Matrix 
Procedure by Bartlett and Kendall (Ref. 16) provided a major 
breakthrough in the field of fast computer solution of the 
steady flow boundary layer problem. The method was developed 
as a means of obtaining rapid computations in the laminar 
regime for steady, multicomponent, chemically reacting 
boundary layers; and the logical extension of the technique 
was to encompass the unsteady flow condition. 

The purpose of this investigation is to develop a com- 
puter program capable of solving the unsteady, laminar 
boundary layer problem by adapting the rapid solution tech- 
niques of the Boundary Layer Integral Matrix Procedure . 
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II. ANALYSIS 



A. GOVERNING EQUATIONS 

In this development, utilizing the integral matrix 
method, it is convenient to make the assumption of incom- 
pressible isothermal flow. The resulting laminar boundary 
layer equations for two-dimensional (n=0) or axisymmetric 
(n=l) flow are as follows: 

Continuity : 



n 9s 3y 



( 1 ) 



Momentum: 



9u ^ ^ = _ i 2? + Ji 

3t ^3s ^9y p 3s p 3^2 



( 2 ) 



,,3U _ 1 3P 

3t '^3t “ p 3s 



(3) 



where s and y are respectively the streamwise and normal 
components, u and v are the corresponding velocity compo- 
nents, U is the streamwise velocity component outside the 
boundary layer, r^ is the radius of the body perpendicular 
to the axis of revolution for an axisymmetric body, n is 
zero for a two-dimensional body and unity for an axisymmetric 
body, t is time, p is density, y is viscosity and P is 
pressure . 
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B. BOUNDARY CONDITIONS 



The boundary conditions considered for the governing 
equations as stated are: 



u 



w 



0 

0 




U(s,t) 

0 



(4) 

(5) 

( 6 ) 

(7) 



where the subscripts w and e refer respectively to wall and 

boundary layer edge conditions. Equation (4) is the require- 
ment of zero slip at the wall, equation (5) implies no mass 
transfer at the wall, and equations (6) and (7) are repre- 
sentative of the assumption that the boundary layer edge 
velocity is equal to the free stream velocity. 

The free stream velocity may be written: 



U = (U^ - b^s) (1 + e cos tot) 



( 8 ) 



where is the mean free stream velocity at the leading 
edge, b^^ is the rate of change of the mean free stream 
velocity with distance from the leading edge, e is the ratio 
of the amplitude of the free stream oscillation to the local 
mean free stream velocity and to is the frequency of the 
oscillation . 
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C. TRANSFORMATION OF EQUATIONS 

Paralleling the development of Bartlett and Kendall 
(Ref. 16) a transformation of variables is performed util- 
izing a combination of the Levy and Mangier and the Howarth- 
Dorodnitsyn transformations which is known by a variety of 
names including Levy-Lees , Lees-Dorodnitsyn and Mangler- 
Dorodnitsyn. This transformation is as follows: 



o 

r ^u y 

n = ^ ^ / p dy (10) 

/2T o 

Simplifying equations (9) and (10) to encompass the 
assumptions of incompressibility and constant temperature 
the transformations become: 






s 

py f 
o 



u r ds 
e o 



( 11 ) 



n = 




( 12 ) 



The transformation to Levy-Lees variables is greatly 
simplified by use of the operator: 






where f is the stream function defined as: 
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( 14 ) 



f - f = / — dn 

w u 

o e 



the prime represents partial differentiation with respect to 
ri so that: 



f' = i 

u 

e 



(15) 



and ot* is a normalizing parameter: 



puu r 



e o 






(16) 



In addition, the partial derivatives are: 



= a*r /2C {iM- + It- 



9n 



(17) 



9( ) _ ot* , . 
9y P ^ ^ 



(18) 



O" < )' 



(19) 



Applying equation (15) to the velocity in the partial 
derivative with respect to time yields: 



3u _ 8f e 

3t “e 3t 3t 



( 20 ) 



Using equations (13) through (20) and applying equation 
(6) to equation (3) results in the following transformed 
momentum equation: 
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Ill i( 1 2 

f + ff + 3(1 - f ) 



-2 {f 



9f 

ding 



-f 



af 



}- 



py 

..2 






3£nu 

e 

at 



-}=0 ( 21 ) 



where 



3 



a£nu 

2 — 

ding 



( 22 ) 



Transforming the boundary conditions, 
through (7) become: 



f 

f 



w 



w 



0 

0 



f 

e 



f 

e 



1 

0 



equations 



(4) 



(23) 

(24) 

(25) 

(26) 



D. TAYLOR SERIES EXPANSIONS 

At any streamwise station s, the boundary layer may be 
divided into N-1 strips which are joined by N nodal points 
where i is 1 at the wall and N at the boundary layer 
edge. The primary dependent variables and their normal 
derivatives are related by Taylor series expansions such 

I 

that the are represented by connected cubics with con- 
tinuous first and second derivatives at the nodal points. 
This is commonly called a spline fit. 

Expanding the primary variables in Taylor series and 

I II 

truncating at f^ the following linear equations result: 
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-f.,, + f.+f: 6n H- f'.' o 

1+1 1 1 ' 1 2 18 1+1 24 


(27) 


1 1 11 III . 2 

-fi +1 + fi + fi £n + fi ^ ^ = 0 


(28) 


II 11 III III 

-fi+i " fi + ^ "i+i T = ° 


(29) 


6n = rii+i - Tii 


(30) 


III III 

i i+1 6n 


(31) 



Equations (27) through (29) represent a set of 3(N-1) 

I 

equations where f is represented as a quartic, f as a cubic^ 
f" as a quadratic and f'" is considered to vary linearly 
between T]^ and 

E. AXIAL AND TIME DERIVATIVES 

The axial derivatives of equation (21) are treated as 
logarithmic two or three-point differences such that: 



2 \ILI\ = ^ h-2 

where for two-point differences: 


(32) 


d = — d, = -d d. = 0 

o 1 o 2 


(33) 
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and for three-point differences: 



<^0=2 A , A„ , '^1 = "2 „A„ , „ ,A 



a 1-2 






"^2 ^. A „ 



u-1 



11-2 H -1 1-2 
(34) 

In equations (32) through (34), ( ) j^^^implies the previous 

streamwise station. 






5,-1 






(35) 



Differencing the time derivatives in the same way: 



l ~3( ) 

L 



= ^o< >k + >k-l ^2( >k-2 



(36) 



where for two-point differences; 



T = -r 
o ,.A 



k“k-l 



Tt = -T 
1 o 



T 2 = 0 



(37) 



and for three-point differences: 



^ k\-lW-2 

o , A, , , A, ^ 1 , A. 



kV2 



k\-l 



k“k-l k“k-2 "• k“k-l k-l*k-2 ^ k\-2 k-l'^k-2 



(38) 



In equations (36) through (38) , ( )j^jimplies the previous 

time . 



k\-l ^k ■ “^k-l 



(39) 
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Equations (37) through (39) may be greatly simplified 
by considering the utilization of equal time spacings . The 
resulting equations for the treatment of time derivatives 
over equal time spacings are as follows: 

Two point differences: 



T 



o 



6t 



h = 



-T 

o 



= 0 



(40) 



Three-point differences: 




1 

26t 



(41) 



where : 



- ^k-i = ^k-i 



^k-2 



(42) 



F. INTEGRATION OF THE BOUNDARY LAYER EQUATION 

The transformed streamwise momentum equation, equation 
(21), is integrated at constant C from to r\^, yielding 

r ..1i i M i ,2 

f + / ff dn + 3 / (1 - f ) dn 

L Ji-1 i-1 i-1 



1 

-2 / 

i-1 



(f 



8f 

ding 



- f 



9f 

ding 



•)dn 



py r J 

*2 at 

a* 1-1 



- ^ ; (f 

a* i=l 



dinu 

1) dn = 0 



(43) 



Expanding equation (43) in Taylor series and applying 
equation (36) to the time derivative terms: 
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(44) 



, 3f ' _ i 8f 

at I at 



i-l 



^ ^ 1^-1 ^ ^ 2 ^ k - 2 ^ 



i ailnu ailnu ^ 

^ <^'- -8t^ = -Tt^ ^L^J. - 



ailnu 



i-l 



i-l 



= {["fl -6n) [T ilnu + T, ilnu +T^ilnu 1(45) 

-'i-l ° ®k 1 ®k-l 2 e^_2 



the integrated boundary layer equation becomes: 






-f[T ilnu + T, £nu + £nu ] 
o 1 e,^_^ 2 e ^_2 



1 

i-l 



■(1 + 3 + 2d^) [f^ XM^ + XM^ + XM3 + 



-2[f. ZM, + f. ZM^ + f. ZM, + f. .ZMJ 
'■ 1 1 1 2 1 3 1-1 4 



+ 6n{3 + [T ilnu + T,5,nu + T^ilnu ] } = 0 (46) 

a*^ ° ®k ®k-l ^ ®k-2 



The Taylor series expansions, XM and ZM, of the integrals 
in equation (43) are listed in Appendix A. 

G. RECURRENCE FORMULAS 

The boundary layer equations, equation (46), and the 
Taylor series expansions of the primary variables, equations 
(27) through (29), together with the boundary conditions. 
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equations (23) through (26) , are solved by Newton-Raphson 
iteration of 4N simultaneous equations in 4N primary 
variable unknowns. These are summarized as follows: 



Primary Variables 



No. of Equations 



No. 



of 



f 

f,', 

Eqns . 



Req : 



N 

N 

N 

N 

4N 



Taylor Series Expansions 
Boundary Layer Equations 
Boundary Conditions 
No. of Eqns. Available: 



Equation Nos . 



No. of Equations 



(27)-(29) 

(46) 

(34)-(26) 



3(N-1) 

N-1 

4 

4N 



The Newton-Raphson solution technique is illustrated by 
considering two simultaneous nonlinear equations: 



F(x,y) = 0 



G(x,y) = 0 



(47) 



which have the solution x=X and y=Y. Defining x and y as 

m ‘‘m 

the values of x and y for the m^^ iteration, the desired 
solution f(X,Y) is expressed in the Taylor series expansions 

3F(x„,y„) 



0 = F(X,Y) = F(x„,y^) + (X - x„) 
m m m 



9x 



8F(x„,y„) 

+ 3y ^ 



0 = G(X,Y) = G(x„,y^) + (X - x„) 
m m m 



9G(x^,y^) 
m m 

9x 



3G(x ,y ) 

+ - V 3y ” ^ ••• 



(48) 
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In the Newton-Raphson method (X,Y) is replaced by 

) and the nonlinear terms in x and y appearing in 
equations (48) are ignored yielding the set of recurrence 
formulas : 



m 






Ay 

■^m 


9F(x ,y ) 
m'-^m 




9x 




9y 


Ax 

m 


8G(x ,y ) 
m m 




m 


9G(x ,y ) 
m'^m 




9x 




ay 


where : 












m 


^m+1 






= 

m 


^m+l - 


The Ax 

m 


and Ay are 
■^m 


corrections to 


be added 



(49) 



(50) 



m m 

respectively yielding and 

are the values of the original functions F(x,y) and G(x,y) 

evaluated for the value of the variables during the m^^ 

iteration. Noting that F (x ,y ) and G(x ,y ) approach zero 

m^m m^'m 

as the corrections approach zero, it is appropriate to con- 
sider them as errors associated with the original equations 
(47) . 

Differentiating the Taylor series expansions, equations 
(27) through (29) , with respect to the primary variables in 
accordance with the derivation of equation (49) yields the 
following recurrence formulas for the m^^ iteration: 

(-1 ) Af^^^t (1) Af^+ (6n ) Af (^) Af% (^) Af ^'+(^) Af^" ^=- ERROR 

(51) 
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(-1) Af (1) Af]^+ (6n) Af % (^) Af^" + 






ERROR (52) 



(-l)Af^^^+(l)Af" + (^)Af" + (^)Afi+i =- ERROR (53) 



where Af^ represents the correction for the coefficients 

of the corrections represent the partial derivatives of the 
Taylor series expansions with respect to the primary var- 
iables, and the ERRORS are equations (29) through (31) 
respectively, evaluated for the values of the variables 
during the iteration. 

The recurrence formulas for the boundary layer equations 

are : 



j^Af" + Af'(l+d^)f + f'(l+d^)Af + Af'(d^fj^_^+ 




1i 

{T Af - Af (T £nu + T,£nu + T^£nu )}| 

° ° ^ ®k-l 2 e^_2 



-2(1+ 3 + 2d^) {Af^XM^+Af ^XM2 + Af^" XM^+ Af^”j_XM^} 
-2{Af^ZMj^ + Af^ZM2 + Af^" ZM^ + Af^_^ZM^}= - ERROR (54) 



where ERROR is the left-hand-side of equation (46) evaluated 
for the values of the variables during the iteration. 

Finally, the recurrence formulas for the boundary con- 
ditions, equations (23) through (26), are: 



Af' = 
w 


- ERROR = - 


(f ' ) 
w m 


(55) 


Af = 
w 


- ERROR = - 


w m 


(56) 
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( 57 ) 



Af ' = - ERROR = 1 - (f ' ) 
e e m 

Af^ = - ERROR = - (f^)^ (58) 

Equation (54) may be further refined by collecting terms 
and writing the equation in terms of correction coefficients 
as follows: 

ClAf^ + C2Af^_^ + C3Af^ + C4Af^_^ 

+ C5AfV + C6Afy_^ + C7Af”' + C8Af^"^= - ERROR (59) 



The correction coefficients Cl through C8 are given in 
Appendix B, 



H. MATRIX FORMULATION 

The coefficients for the recurrence formulas, equations 
(51) through (53) and (55) through (59) , evaluated for the 
m^^ iteration form a square matrix [^aJ with 4N rows and 
columns. The matrix equation is: 



A 


I 

> 

< 

1 


i 


['] 


_ 


L 







where the column matrix LAVj is the matrix of corrections to 
the primary variables, and the column matrix |*eJ is the 
matrix of errors. The matrix ^aJ contains many zeros, and 
great saving^ in computer storage space and computation time 
may be realized through proper ordering of the equations and 
variables. 
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The equations are first divided into linear and nonlinear 
sets, and at the same time, the variables are classified as 
linear or nonlinear, resulting in the partitioned matrix 
equation : 



AL 


BL 




AVL 


= - 


EL 


ANL 


BNL 




AVNL 


ENL 



The division of variables into linear and nonlinear sets is 
generally arbitrary except that the division must be made to 
insure that the square matrix AL is not singular. Adopt- 
ing the variable groupings of Bartlett and Kendall (Ref. 16) 
the linear and nonlinear corrections are arranged as follows: 



AVL (Af2 ,Af^, . . . ,Af^,Af^,Af2 , . . . ,Af^,Af^,Af^ , . . . 



AVNL (Af ,Af " , Af ") Af . . . , Af '" ) 
w' e' w' 2' e 



The linear and nonlinear equations are sequenced within the 



partitioned square matrix j^Aj as shown in Figures 1 and 2, 
which demonstrate the case of 4 nodes across the boundary 
layer. 

Expanding equation (61) , the resulting linear and non- 
linear equations are, respectively: 



AL J ^ AVL j + I BL [ 



1 

avnlJ = 



-EL 



r ANL 1 


1 avl] 


+ BNL 


: AVNL 1 = -ENlI 


L J 


1 L J 


! L 


1 , J w J 



(62) 

( 63 ) 



Solving equation (62) for the linear corrections: 
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FIGURE [.LINEAR PORTION OF [A] MATRIX 
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r -| r nr 1-1 ~ 

AVL = -' AL i BL ; 1 AVNL i + ^ AL ! ; -EL 



j L 



( 64 ) 



r 1 r 

\ p- 



Examination of the j AL j and '■ BL ; matrices reveals that 

-1 - ^ -i 

the product f AL~j T BL~j may be performed given only the nodal 
spacing and must be determined only once for a given problem, 



thus : 



where : 



j AVL 


- - 


\ 





BAlJI AVNL 






ELA 



-» •' -j 

, BAl ' = i AL I 

I . ' I 



-If -] 



ELA 



I AL 

^ J 






(65) 

( 66 ) 

(67) 



Applying equation (65) to equation (63) , and solving for 
the nonlinear corrections: 

1 ^ 1-ir' 



AVNL : = : BNL 



ENL 



where : 



^ n - 

BNL - = BNL ; - ANL 

■ ■' ' "^1 r " 

ENL = : -ENL J - j ANL 



BAl 



I^elaJ 



( 68 ) 

(69) 

(70) 



I. SOLUTION PROCEDURE 

The matrix solution proceeds according to the following 
algorithm: 



1. Given the nodal spacing determine r-^Ll 
the product j^BAl~|. 



and form 



2. Given initial values of the primary variables deter- 
mine -EL^ and T-ENlJ, and form the product j^ELAj. 
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r -} f r 

3. Form. I BNL , and , ENLj from the coefficients of | ANL 
and r BNL \ ‘ 

-Ir- 



4. Invert TbnlI and form the product 
which is the solution for rAVNL"|. 



BNL 



r 



ENL 



5. Solve for UvLj from equation (64) using I AVNLJ 
as determined in srep 4 . 



6. Add the linear 






to the corresponding primary 



variables to complete the iteration. 

7. Compute the new errors T-ElJ and 
test for a maximum allowable error, if this criterion 
has not been reached, return to step 2 using these 
errors to commence the new iteration. 



T-enl], 



and 
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III. RESULTS AND DISCUSSION 



The computer program has been used to investigate three 
basic categories of boundary layer problems using a flat 
plate model; steady Blasius flow (zero pressure gradient), 
steady Howarth flow (uniform adverse pressure gradient) , 
and oscillating Blasius flow. The velocity profiles which 
resulted from these solutions are compared to previously 
reported results and presented in Figures 3, 4, and 5. 

A. STEADY BLASIUS FLOW 

In Figure 3, computed values corresponding to three 
nodal spacing choices are superimposed on the Blasius pro- 
file. Agreement between the numerical solution and the 
classical solution is better than 0.08 percent for as few 
as 10 nodal points. 

B. STEADY HOWARTH FLOW 

Results for a steady flow in a uniform adverse pressure 
gradient are compared with the solution of Howarth (Ref. 11) 
in Figure 4. The comparison of the computer solution with 
that of Howarth is accomplished by correlating the dimen- 
sionless parameter x* (=b^s/U^) , which characterizes 
Howarth 's velocity profiles. Agreement between the two 
solutions is generally good with differences no greater than 
4.35 percent for 22 nodal points. The differences encoun- 
tered may be reduced somewhat by increasing the free stream 
velocity of the test problem. This results in larger 
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VELOCITY PATIO {f = u/ug) 

FIGURE 3. STEADY BLASIUS FLOW DATA COMPARISON 
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6.0 
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VELOCITY RATIO (f'=U/Ug) 

FIGURE 4. STEADY HOWARTH FLOW DATA COMPARISON 



stream-wise pressure differences in the equations of motion 
and consequently lesser roundoff errors. 

C. OSCILLATING BLASIUS FLOW 

The unsteady capabilities of the numerical procedure 
were investigated in the solution of Blasius flow with 
sinusoidal oscillations superimposed on the mean flow. 
Utilizing the critical oscillation frequency (=.6U^/s, 
for a flat plate) of Lighthill's analysis (Ref. 2) to 
determine the solution approximation with which to compare 
the numerical results, the profile of velocity fluctuations 
which resulted was compared with the "quasi-steady" 
analysis as shown in Figure 5. Agreement between the calcu- 
lated velocities and the Lighthill profile is within 15 
percent deep in the boundary layer and within 5.5 percent 
in the outer regions. Noting that Lighthill's analysis is 
based on linearized small perturbation theory and that the 
amplitude of oscillation in the test problem of 0.1 can- 
not be considered a small perturbation, better agreement 
could scarcely be expected. 

D. NUMERICAL PROCEDURE 

The computer program has been operated on both IBM 360/67 
and CDC 6600 computer installations. Core storage require- 
ments are 130,000 bytes for a maximum array size encompassing 
25 nodal points across the boundary layer, 50 streamwise 
stations and 50 time increments per oscillation cycle. A 
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RATIO OF VELOCITY FLUCTUATIONS (AUj/aU,) 



FIGURE 5. OSCILLATING BLASIUS FLOW 
DATA COMPARISON 
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typical 25 nodal point Blasius solution was performed in 
7 seconds on the CDC 6600. 

Limitations in the present numerical capabilities of 
the program precluded a more extensive investigation of 
parametric effects, however the results qualitatively 
demonstrated the feasibility of the computational procedure. 
Further refinements to the computer program should enable 
a more extensive investigation to include quantitative com- 
parisons with experimental results. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



The feasibility of the Unsteady Boundary Layer Integral 
Matrix Procedure in the treatment of simple steady and 
unsteady flow problems has been demonstrated. Favorable 
comparison with classical results provides a firm base for 
continued investigation in the unsteady regime. 

Continued research utilizing the program should include 
the following revisions: 

1. Alteration of the present time derivative component 
arrays (FTMl) and FTM2) to insure that all values of the 
components are stored for each nodal point, instead of 
retaining only the values at the wall. 

2. Incorporation of Blasius profiles as first estimates 
for the primary variables in order to decrease the 

time required to complete the first nodal iteration 
procedure . 

3. A study of array usage aimed at reducing the re- 
quired storage and subsequent conversion to double 
precision computation in order to minimize round-off 
interference . 

4. Addition of the capability of handling compressible 
flows in order that higher velocities may be introduced. 

5. Inclusion of turbulent boundary layer capabilities 
perhaps through the use of eddy transport properties. 
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APPENDIX A 



TAYLOR SERIES EXPANSIONS OF INTEGRALS (Ref. 16) 

i r ^ 9 

/ ff" dri = I ff ' j - / f ' dn 

i-1 L Ji-1 i-1 



I 9 I II III IM 

/ f dn = + ^i ™2 ™3 ^i-l ™4 



»'i - »" < '1 - -t" - ‘I" ^ * 'r:i #1 



XM, .-Sn^(l f! - tr t t"', 

2 ''2 1 13 1 120 1-1 30 



„„ _ ^"U6n , .... W , 56n^ 

XM3 - 6n (g f . - f ^ 1^ + ^i-i ToT 



XM, = 6n^^f: - fV + f"' W ^ f". 

4 24 1 1 30 1 504 i“l 252 



1 

2 / (f 

i-1 



• df _^i. 3f 






)dn = -[d^ff' + + d 2 f,_ 2 f 



+ 2d [f!XM^ + fVxM. + f!”XM- + f."' XMJ 
Oil 1 2 1 3 1-1 4 

+ 2[f!ZM-, + + f!"zM- + fl",] 

'■ 1 1 1 2 1 3 1-1 



2 2 

ZM^ = 6 n ( YM^ - YM 2 -^ + YM 3 ^ + YM^ 



ZM 2 =-6n2 (i YM^ - YM 2 ^ + YM 3 + 



6n" 
"4 30 



YM, 

2 r 2 



ZM3 = 6n^| ym^ - ym2^ + YM3 ^ + YM4 



(Al) 

(A2) 

(A3) 

(A4) 

(A5) 



(A7) 

(A 8 ) 

(A9) 

(AlO) 
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(All) 



,3 , 1 



6n 



ZM3 = 6n (^YM^ - YM2 ^ + YM3 + ...4 — 



W . YM. ^ 



™1 = <^ 1 ^ 1,1 ^ '^ 2 ^ 2,1 



™2 “ ‘^l^il-l,i ‘^2^«.-2,i 



YM, = d,f"', . + d,f. , . 

3 1 i-1,2. 2 £-2,1 



™4 ^l^£-l,i-l ^2^£-2,i-l 



(A12) 

(13) 



(A15) 
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APPENDIX B 



BOUNDARY LAYER EQUATION CORRECTION COEFFICIENTS 



Cl = (1+d )f!-£^ {T (l-dnu )-(T, )lnu + T,£nu )} 
o 1 o e 1 ^ %-2 



C2 =-(l+d )f.' + -Sir {T (1-S.nu )-(T,£nu + T.S,nu )} 

o 1-1 ° el 2 



C3 = (l+d^)f.+(d^f,_i^i+d2f,_2,i)-2(l+6+2d^)XM^-2ZM^ 






C5 = 1 - 2(l+B+2d )XM- - 2ZM- 
O 2 2 



C6 = -1 



C7 = -2 (l+3+2d^)XM^ - 2ZM^ 



C8 = -2(l+3+2d )XM. - 2ZM. 

O 4 4 



(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 

(B8) 
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APPENDIX C 



PROGRAM DESCRIPTION 



A. GENERAL 

The computer program referred to as UBLIMP, is written 
in FORTRAN IV source language and has been operated on both 
IBM 360/67 and CDC 6600 computer installations. Core 
storage requirements are 130,000 bytes for a maximum array 
size encompassing 25 nodal points across the boundary layer, 
50 streamwise stations and 50 time increments per oscillation 
cycle . 

The program consists of a main program (UBLIMP) and 21 
subroutines which are divided into two sections as shown 
in Figure 6. LINK 0, common to both sections, contains the 
main program and certain service routines which are common 
to both sections in the fully expanded program. LINK 1 
sets up the boundary conditions and LINK 2 iterates to solve 
the boundary layer equations. 

The program may be applied to steady or unsteady flow, 
two-dimensional cartesian or axisymmetric bodies with similar 
or nonsimilar profiles. It allows quadratic or cubic curve 
fits of the primary variables. Program options are con- 
trolled by the control variable KR and are detailed in the 
input instructions for the UBLIMP program in Appendix D. 

A program linkage schematic is shown in Figure 7 and a 
functional flow chart is shown in Figure 8. Correlation 
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— UBLIMP 



LINK 



LINK 0 



— ZEROIT 

— MATONE 

— SLOPQ 

— TAYLOR 

— RERAY 



— SETUP 
— RECASE 
— INPUT 
— STDETA 

1 — LINMAT LINK 2 

— FIRSTG 
— REFCON 
— HISTXI 



FIGURE 6. COMPUTER PROGRAM OVERLAY 



— ITERAT 
— LINGER 

— ABMAX 

— MONGER 
IMONE 

— AMSET 

— lONLY 

— OUTPUT 

STRUCTURE 
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FIGURE 7. SCHEMATIC OF UBLIMP SUBROUTINES 



CASE 




FIGURE 8. FUNCTIONAL FLOW CHART OF UBLIMP PROGRAM 
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between the two is indicated by letters running from A 
through K. 

B. PROGRAM FUNCTIONS 

1 . Main Program (UBLIMP) 

The main program provides linkage between the 
boundary conditions section and iterative solution section 
of the program. Initial core zeroing and the call for 
solution output are also initiated by the routine. 

2 . Subroutine ZEROIT 

ZEROIT performs the initial zeroing of all named 
common blocks to eliminate the possibility of introducing 
random numbers in the solution. 

3. Subroutine SETUP 

SETUP is the control routine for setting up the 
boundary layer edge conditions and the streamwise derivatives 
for a new time, station or case. 

4 . Subroutine RECASE 

RECASE is the control routine for the input of 
boundary layer data and initializes control variables for 
the treatment of surface discontinuities. 

5. Subroutine INPUT 

INPUT reads all the data for a problem into the 
computer. 

6 . Subroutine STDETA 

STDETA permits simplification of the input data deck 
by providing standard nodal point distributions should the 
user not care to provide one. 
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7 . Subroutine LINMAT 

LINMAT sets up the matrices for the Taylor series 
expansions and linear boundary conditions from the eta 
spacing and solves to express the linear corrections in 
terms of the nonlinear corrections. 

8 . Subroutine MATQNE 

MATONE performs operations on a column of a matrix 



of coefficients or on a column of a matrix of errors 



^-ElJ (designated E in the subroutine call list) such as to 
form [^j* 

9 . Subroutine FIRSTG 

FIRSTG computes first estimates of the primary 
variables based on the nodal point distribution should 
the user not care to provide them. 

10 . Subroutine REFCON 

REFCON calculates the boundary layer edge conditions 
and sets up the wall boundary conditions. 

11 . Subroutine SLQPQ 

SLOPQ defines a cubic equation for a set of points, 
calculates the average slope at each point and integrates 
the cubic equation between each pair of points. 

12. Subroutine HISTXI 

HISTXI computes terms involving the axial deriva- 
tives and time derivatives and stores those upstream quan- 
tities needed for these difference relations. 
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13. Subroutine TAYLOR 

TAYLOR calculates the coefficients in the Taylor 
series expansions of the integrals which appear in the 
integral form of the boundary layer equations. 

14 . Subroutine ITERAT 

ITERAT is the control program for performing the 
boundary layer iteration. 

15. Subroutine LINGER 

LINGER evaluates errors for the linear equations 
and determines the maximum linear error. 

16 . Subroutine ABMAX 

ABMAX searches a given array for the entry with 
the maximum absolute value. 

17. Subroutine NONCER 

NONCER performs that portion of a boundary layer 
iteration having to do with the solution of the nonlinear 
corrections of the nonlinear equations, computes the damp- 
ing factor and applies it to the corrections and then 
applies the damped corrections to the primary variables. 

18 . Subroutine IMONE 

IMONE evaluates the coefficients of (i-1) cor- 
rections for the i^^ nonlinear equations, where i is the 
i^^ nodal point in the boundary layer. 

19 . Subroutine AMSET 



AMSET calculates the contributions of the ANL and 



applies them 
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20. Subroutine lONLY 



lONLY evaluates the coefficients of the i^^ cor- 
rections for the i^^ nonlinear equations, where i is the 
i^^ nodal point in the boundary layer. 

21 . Subroutine RERAY 

RERA.Y replaces a matrix fcj with its inverse and 

r 

forms the product of the inverse with another matrix | D .. 

22. Subroutine OUTPUT 

OUTPUT prints a standard boundary layer output block 
for a converged solution, or at the end of each iteration 
if desired. 
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APPENDIX D 



INPUT INSTRUCTIONS FOR THE UBLIMP PROGRAM 



GROUP 1 CONTROL CARD, TITLE 
CARD 1 FORMAT(20I1,15A4) 

FIELD 1 (Columns 1-20) This is the variable KR which is 
used to control the various program options. 

COLUMN 1 Determines whether a new set of ETA values 
is to be input for the present case. 

0 Uses resident values from previous case. 

1 Values input by user (mandatory for first case). 

COLUMN 2 Designates the type of first guesses to be 
utilized for primary variables. 

0 Uses built-in relations to calculate first 
guesses . 

1 First guesses input by user. 

2 Uses resident values from previous case (cannot 
be used for first case) . 

COLUMN 3 Determines the treatment of streamwise 
derivatives . 

0 Performs similar solution at each streamwise 
station . 

1 Considers two-point difference relations at all 
stations except that a similar solution is per- 
formed at the first station for non-blunt bodies 
and at the first two stations for blunt bodies. 

2 Considers three-point difference relations at all 
stations except that a similar solution is per- 
formed at the first station and a two-point solu- 
tion is performed at the second station for 
non-blunt bodies; similar solutions are performed 
at the first and second stations and a two-point 
solution is performed at the third station for 
blunt bodies; and a two-point solution is per- 
formed for the first station after a discontinuity. 
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COLUMN 4 Determines when output block is to be 
printed . 

0 Output block printed for converged solution or 
for nonconverged solution after 50 iterations. 

1 Output block printed after each iteration. 

COLUMN 5 Determines treatment of Entropy Layer. 

(Not used in present program. KR(5)=4 
is input to maintain input data sequencing.) 

COLUMN 6 Designates body shape. 

0 Axisymmetric blunt body. 

1 Planar blunt body 

2 Axisymmetric sharp body. 

3 Planar sharp body. 

4 Axisymmetric or planar shape which has no sharp 
tip or blunt stagnation point, such as a nozzle. 

COLUMN 7 Steady, zero pressure-gradient noise filter. 

0 Normal computations of DUDS. 

1 DUDS is forced to zero. 

COLUMN 8 Designates form in which wall mass fluxes 
are input. KR(8) is not utilized if wall 
mass fluxes are not input. 

0 Wall fluxes input in LBS/SEC FT**2 

1 Wall fluxes input in normalized form (divided by 
-ALPHASTAR) . 

COLUMN 9 Together with KR(ll) this designates the 
type of wall boundary conditions. 

0 Assigned stream function at the wall. 

1 Assigned mass flux at the wall. KR(9)=1 is used 
if zero mass flux is used. 

COLUMN 10 Determines the type of curve fits employed 
to represent the primary variables. 

(KR(10)=0 is recommended for accuracy for 
most problems, however for severe problems 
KR(10)=1 is better since the cubics can 
become poorly behaved.) 
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0 utilizes connected cubics. 



1 Utilizes connected quadratics except for the 
outermost segment where connected cubics are 
used. 

COLUMN 11 Together with KR(9) this designates the 

type of wall boundary conditions. (Presently 
KR(11)=0 is the only active value and 
requires an assignment of wall temperature.) 

COLUMN 12 Provides access to standard nodal-point 
distributions built into the program. 

0 User inputs data for nodal-points. 

1 8 points with more concentration at the wall. 

2 8 points equally spaced. 

3 10 points equally spaced. 

4 12 points equally spaced. 

5 14 points equally spaced. 

6 15 points with more concentration at the wall. 

7 18 points equally spaced. 

8 22 points equally spaced. 

9 25 points equally spaced. 

COLUMN 13 Permits the assignment of a convergence 
damping factor. (This factor is over- 
ridden if a smaller damping factor is 
computed internally by some constraint.) 

0 No damping factor is assigned. 

J If J is greater than zero, corrections are 
damped uniformly by J/10. 

COLUMN 14 Non-zero entry causes a complete set of 

primary variables to be output for future 
use as first guess inputs. 

COLUMN 15 Non-zero entry provides debug output for 
first guesses and linear matrices. 
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0 No debug. 

1 First guesses are dumped. 

2 Linear matrices before and after inversion are 
also dumped. 

COLUMN 16 Non- zero input indicates the pressure pro- 
file has been input in the form of pressure 
coefficients . 

COLUMN 17 Non- zero entry provides debug output for 
coefficients . 

0 No debug. 

Y For (Y+l-ITS) greater than zero, where ITS is the 
number of the current boundary layer iteration, 
the coefficients which combine to make up the 
nonlinear equations are dumped and the derivatives 
of the nonlinear equations with respect to the 
nonlinear variables are dumped before and after 
inversion . 

COLUMN 18 Not used. 

COLUMN 19 Non- zero entry provides debug output for 
linear and nonlinear errors. 

COLUMN 20 Determines which of the two internal sets 

of first guesses are to be used if the user 
does not provide them. 

0 Straight line Blasius inputs are calculated for 
the stream function and velocity ratio, and zero 
is input for the first and second derivatives of 
the velocity ratio. 

1 Straight line Blasius inputs are calculated for 
all the primary variables. 

FIELD 2 (Columns 21-80) CASE. Title of the case 

(ALPHA-NUMERIC) . Used for identification of 
printed output. 

GROUP 2 TIMES AND STATIONS 
CARD 1 FORMAT (13) 

FIELD 1 (Columns 1-3, right- j ustified) NITEM Number of 
time grid points per cycle when considering an 
unsteady solution, otherwise input 1. (max=50) 
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CARD 2 


FORMAT (ElO. 4) 


FIELD 


1 (Columns 1-10) TIME(l) Value to be used in 

the output block to identify a time of solution. 


CARD 3 


FORMAT (13) 


FIELD 


1 (Columns 1-3, right- justified) NS. Number of 
streamwise stations. (max=50) 


CARD 4 


FORMAT (8E10. 4) 


FIELD 


1,2,... 8 per card (10 columns per field) S(IS) 
Streamwise distance upon which the boundary 
layer solution is based in feet. Blunt body 
problems should start with S=0 . 0 . Sharp body 
or nozzle problems must start with some finite 
value of S. The boundary layer is assumed 
similar up to and including this first station. 
A negative entry for S indicates a discon- 
tinuity at that station, and alters the differ- 
encing scheme for axial derivatives. 



GROUP 3 NODAL DATA (Skip this GROUP for KR(12 )t>^ 0) 



CARD 1 


FORMAT (13) 


FIELD 


1 (Columns 1-3, right- justified) NETA Number of 
nodal points across the boundary layer, 
(maximum of 25) 


CARD 2 


FORMAT (8E10. 4) 


FIELD 


1,2,... 8 per card(10 columns per field) ETA(I) 

ETA stations across the boundary layer starting 
at the wall (ETA=0.0). It is recommended that 
a value between 5.0 and 6.0 be assigned at the 
boundary layer edge. UBLIMP values of ETA are 
equivalent to Blasius values divided by the 
square root of two. 



GROUP 4 BODY SHAPE DATA 

CARD 1 FORMAT (2E10. 4) USED ONLY IF KR(6)=0 or 1. 

FIELD 1 (Columns 1-10) CONE. Cone half angle in sphere- 
cone shape bodies. Leave blank for other body 
shapes . 

FIELD 2 (Columns 11-20) RNOSE. Effective nose radius 
in feet. Used to calculate stagnation point 
velocity gradient from Newtonian relations. 
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CARD 2 FORMAT (8E10. 4) USED ONLY IF KR(6)=0, 2 or 4 

FIELD 1,2,... 8 per card (10 columns per field) ROKAP(IS) 





This is the local body radius in feet normal 
to the body centerline raised to the n power, 
where n is unity for axisymmetric bodies and 
zero for planar bodies. Therefore, ROKAP is 
unity for planar bodies and the local body 
radius for axisymmetric bodies. For planar 
bodies, this card set is used only if KR(6)=4. 

(A special input format can be used for spher- 
ical-nosed bodies as follows. Set ROKAP (1) 
equal to minus the nose radius. The nose radius 
is then set to -ROKAP (1) and ROKAP (1) is set to 
zero. If subsequent ROKAP are input as zeros, 
the program computes ROKAP from S for a spher- 
ical nose. The first non-zero entry is then 
ROKAP at the tangent point. If this again 
followed by zeros , linear interpolation is used 
to the next non- zero entry to yield ROKAP along 
a conical afterbody.) 


GROUP 5 


STAGNATION DATA 


CARD 1 


FORMAT (ElO .4) 



FIELD 1 (Columns 1-10) PTET. Local stagnation pressure 
in atmospheres. 

GROUP 6 FIRST GUESS INFORMATION 
CARD 1 FORMAT(5 (2XE14.7) ) 

FIELD 1,2,... 5 per card (16 columns per field) F(I,J). 





First guesses for stream function (F(1,J)), 
velocity ratio (F(2,J)), shear function 
(F(3,J)), and derivative of shear function 
(F(4,J)) across the boundary layer. 


GROUP 7 


STREAMWISE DISTRIBUTIONS FOR EDGE CONDITIONS 


CARD 1 


FORMAT (5 (2XE14. 7) ) 



FIELD 1,2,... 5 per card (16 columns per field) PRE(IS) 



CARD 2 


Ratio of local static to stagnation pressure 
for the mean flow condition. When KR(16) is 
non-zero PRE is input in the form of negative 
pressure coefficients referred to PTET. 

FORMAT (Blank) This card is inserted to maintain 
data continuity. In more complex formulations 
this card corresponds to a control flag for 
updating the pressure profile during the com- 
putation . 
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GROUP 8 STREAMWISE DISTRIBUTIONS FOR WALL INPUT CONDITIONS 



CARD 1 FORMAT (ElO. 4) 

FIELD 1 (columns 1-10) T. Wall temperature in degrees 
R (for the present this is the free stream 
static temperature) . 

CARD 2 FORMAT (8E10. 4) USED ONLY FOR KR(9)=0 

FIELD 1,2,... 8 per card (10 columns per field) FW(IS). 
Wall stream function (negative for mass 
addition) . 

CARD 3 FORMAT (Blank) USED ONLY IF CARD 2 IS USED. This 
card is inserted to maintain data continuity. 
In more complex formulations this card cor- 
responds to a control flag for updating the 
wall stream function profile during the compu- 
tation . 



CARD 4 FORMAT (8E10. 4) USED ONLY FOR KR(9)=1 

FIELD 1,2,... 8 per card (10 columns per field) RHOVW(IS). 

Total mass flux at the wall (LB/SEC FT**2 or 
dimensionless for KR(8)=0 or 1 respectively). 
These values are positive for mass injection. 

CARD 5 FORMAT (Blank) USED ONLY IF CARD 4 IS USED. This 
card is inserted to maintain data continuity. 

In more complex formulations this card cor- 
responds to a control flag for updating the 
total mass flux profile during the computation. 

GROUP 9 FREE STREAM DATA 
CARD 1 FORMAT (3E1 0.4) 

FIELD 1 (Columns 1-10) UMFS . Mean free stream velocity 
in FT/SEC. 

FIELD 2 (Columns 11-20) FREQ. Frequency of mean free 
stream oscillations in HZ. 

FIELD 3 (Columns 21-30) EPS. Ratio of the amplitude of 
the free stream oscillation to the local mean 
free stream velocity. 

LAST CARD FORMAT (Al) The purpose of this entry is to per- 
mit a test on whether or not a new case is to 
follow. In the event a case does not converge 
in the alloted number of iterations, any re- 
maining cards for that case are read and then 
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ignored until either a period or a comma is 
encountered in column 1. A comma in column 1 
indicates that another case is to follow, while 
a period in column 1 indicates that this is the 
last card in the input deck. 

The preceding input instructions should permit the user 
to fully exercise the UBLIMP program with a minimum of 
difficulty. For more detailed information on various rou- 
tines the reader is referred to reference 18, which is the 
user manual prepared for the BLIMP program. 
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THETA DELSTAR DELSTAR/ 

(FT) (FT) THETA 

.451E-06 2.190E-05 2.592E 00 
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366E-05 2.400E 00 1.232E 00 9.011E-01 1.676E-01 -2.154E-01 2.703E 01 
822E-03 3.200E 00 1.991E 00 9.801E-01 4.694E-02 -8.629E-02 2.940E 01 
277E-03 4.000E 00 2.784E 00 9.977E-01 6.897E-03 -1.381E-02 2.993E 01 
846E-03 5.000E 00 3.783E 00 l.OOOE 00 -1.935E-05 -2.210E-05 3.000E 01 
416E-03 6.000E 00 4.783E 00 l.OOOE 00 3.293E-08 6.087E-05 3.000E 01 



COMPUTER PROGRAM 



c un«^tfady r"^unda»^y layfr intpgpal matrix procfoure 

r (UPLIMD) 



real MIJ 

CnMMnN/PUMCPM/RUMP,cnRMA ,EASF,ICnRM,WOOT ♦ 1777, IS®, IX 
common /^PGC CM /PPM( 50) ,PTE(50) fODESfU E( 50,50 ), PE (50, 50) 
1,RHD, MU,T,OSIP( 50) ,IDST P 

COMMON' /! NIC OM/KR (20) ,KIN,KntJT,MATlI ,MAT1J,NFTA, I, IS, NS 
1 ,IT,NTIMF,nAM,NLE 0,NNLT0,NRNL,ITS,CASE(15),R(B) ,MWE,NO 
2M, ^TFM,NT TEM,KP17 ,NBT,N8T2 , 1 DENT ,U MFS • FRFO, E®S 
CnMMnN/poMCnM/TiME (50) ,®RE(5Q ) ,PTET,5 ( 50 ) , ROK AP ( 50 ) , RN 
nSF ,VKAP.N'OI SC , ID! cc (50) ,NS0(1G ) ,MSD( 10 ) , I®R=,PAONn, CO 

2ne,delt,(:ycle 

COMMON/PPPCCM/PP ,SC,XM(5) 

CO^^MON/VARCQM/F (5 ,25) ,FTM1 ( 5 , 25 ) ,FT M2 ( 5 , 25 ) 

CALL 7FPCIT 

1 FaP.MAT(Al) 

2 FOPMAKFIO.A) 

3 FOP^AT( 1H1,/////////7/////5X10HRUN TIME = FR.2,RH SECO 

1N0<^ ) 

data lAS^/lH,/ 
data LAST/IH,/ 

KTN=4 

Kni|T=R 



R( 1)=,5 

R( 2 )=. 333333333 
R( 3) = . 166666666 
R(A)=.125 
R( ^)=.CA1666666 
B(^)=,C33333333 
B( 7)=,C1388R888 
B( 8 )=.CC396825A 
RHP=. 07651 
mU=3,71cf-QT 
TT=l 



A 6 MWF=-1 
FA‘^F=1. 

I7^7=C 

IS=1 

41 TTFM=1 
A2 CALL SFTyo 
A3 CALL ItfrAT 
CALL OUTPUT 
IF( MnN)^3,4A,4C 
44 itcm=ITFm+1 



DO 4c T=i,fsjETA 
DO 40 J=l,5 

Ftm 2( J,l )=FTmk j,i) 
FTmK J,!)=F( J,I ) 

4<5 continue 

IF( TTEM-NITFM) 42,42,45 
45 IS=IS+1 

IF( IS-N5) 41 ,41,40 
40 READ(KIN,1) JASt 

TF( lA ST-JAST) 47,46,47 
47 TF(LAST-JAST) 40,43,40 

4P 

END 



SUBROUTINE 7ER0IT 
COMMON/BUMCOM/BUMCn(R ,1 ) 

Common /EDGC0M/EDGC0(5 ,1) ,EDGC('^ ,10 ) ,EDG( 1, 100 ), E0(5, 10 
100) 

CnMMQN/FRRC0M/ERPC0(5 , 1 ) ,E PRC (2 ,10 ) ,ERP(3, 100 ) 

COMMON /ETACOM/E TACO ( <5, 1) ,ETAC(B,10) , ETA ( 2 , 100 ) , ET ( 2, 10 
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ICvO) 

CnMMPN/HT T SCO (6 »1) 

COMMON /!MTCny/i nTC.p( »1 ) 
COMMON /NCNCON/NONC (2 tlO) 
COMMnN/P<^Mr,QW/pp.MCO{ « *1 ) 
CnMMriN'PPPCOM/PRPcn(T,l ) 

COMMON /T<^ ^^C aw/TFMCO( f 1 ) 
COMMON/VAPCOM/VA'^cnC*^ , 1 ) 
COMMON /UiLCOV/WALCO( 3 ,1 ) 
7=0*0 

DO 3 I=1,10CC 
on 1 J=l,2 

FT( j, I )=7 

1 CONTI NUF 
^0 2 J=l,5 
For J,I )=7 

2 CONTI nUF 

3 CONTI NUF 

no 8 1=1,100 

FnO( 1 , I )=7 
WM ( 1 ,t) = 7 

no 4 J=l,2 
J,T) = 7 

PPMf j,i )^7 

^ C^^)TTNIJP 

j=i,3 
FRP ( J, T )=7 
VAP( J,T )=7 

F continue 

DO 6 J=l,4 
HTS( J,U=7 
TFM( J,I )=7 
f CONTI NIT 
DO T j=i ,8 
N0N( J,T )=7 
7 CON'^INIF 
« CONTI NUF 

no 14 1=1,10 

no c j=i,2 

EROC( J,I )=7 
NONC ( J, n=7 
^5 CONTTNIT 
DO 10 J=l,3 
HTSC( J,I )=7 
"^FMc ( j ,n=7 

10 C-ONtimup 

DO 11 J = l,’= 

^nr,C( J,I ) = 7 

11 continue 

no 12 j=l,7 

INTf( .J,I )=7 
nRMC( J,T)=7 
V^PC( J, T )=7 

12 continue 

no 13 J=l,8 

FT^C( J,I )=7 

13 CONTINUE 

14 COMTINUP 

DO 1C j=i 3 
WAlCO( J,lf=7 
IE C'^NtinUP 
DO 20 J=l,6 
HI SCO( J,1 )=7 
20 CONTINUF 
no 16 J=1,E 
pnccri j,n=7 

FRRCO( J,l)=7 
INTCCU,1)=7 
TEMCO( J,I )=7 
VAPCP< J,l)=7 
16 CONTINUE 
DO 21 J=l,7 



) ,HIPC( 3 , 10 ) 
) ,INTC(7,10) 
► ,NnN(8,l00) 
,PRMC(7,10) 

,TFMC(3 ,10) 
,V.^RC(^ , 10 ) 
,WAL( 1 , 100 ) 



,HIS (4, 100 ) 

,PRM(2, 100 ) 

,TFM(4,100) 
,VAP (3, 100 ) 
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PRPCO( J,l)=Z 
21 CONTINUE 
rn 18 J=l,8 
BUMCO( Jtl)=? 
POMCO( J,l)=Z 
18 CONTINUE 
DO IQ J = l,'5 
ET^CO( J,l)=7 
IQ CONTINUE 
RETURN 
END 



SUPPOUTINE SETUP 
REAL MU 

01 MEN SIGN HISTI (278) ,HI ST2 ( 5 1 ) , HIST3 (375) ,HIST4( 100 ),V 
1MAT( 158) 

COMMON /EDGCOM/PEM( 50) .PTE (50) ,DUE$ ,UE(50 ,50 ) , PE ( 50, 50 ) 
1,PH0,MU,T.DSIP(50) , I OSIP 

COMMON /HI SC0M/TZ,T1 ,T2,TF(25) ,C 1 ,C2 , C3, C4, BET A, ZM(5, 25 
1) ,XI( 50) ,HF(25,5) ,HUE ,HHUE , DLX2 ,(^3 M ( 50 ) , BET AM ( 50 ) 
CnMM0N/INTCCM/KR(2 0) , KI N ,K0UT , M ATI I ,MAT1J,NETA, I, IS, NS 
1, IT, NT I ME ,NAM,NLE0,NNLE0,NRNL,TTS,CASE(15 ),B(8), MWE,NO 
2N,ITE»',NITEM,KP17 ,NBT,NBT2 ,I0ENT,UMFS,FRE0,FPS 
COMMON/PPMCOM/TIme (50 ,PRE(50) , PTET,S (50 ) , ROK AP ( 50 ) , PN 
10SE,VK4P,N0ISC,I0ISC(50) ,NS0(10) ,MS0(10), I PR E , RA DNO, CO 
2NE,DELT, CYCLE 

COMMON /VARCCM/F (5,25) ,FTM1 ( 5 ,25 ) , FT M2 ( 5 , 25 ) 
C0MM0N/VIALC0M/FW(50) ,RH0VW(50) , IFW, IRHOVW,.)RHOVW 
EQUIVALENCE (HI STl ,XI ) , (HI ST2,PTE) , (HIST3,»=) , (HIST4, FW 
1) ,( VMAT,TZ) 

2 FORMAT! 1H18X4HTIMEE12. 5, 56H SECCNOS -- - 

I -.6X,A5 ,3X,A2, IH,, A2) 

3 FnRMAT(lHl,3X4HTIME,E12,5,34H$FC0NDS * ♦ ♦ STRFAMWISE 
lOlMENSICN , E12.5,13H FEET 

4 FOPMA -T! 1H1BX4HCA SEE12.5,56H 

X ,6X,A5,3X,A2,1H.,A2) 

5 FORMAT! 1H1,2X12H* * ♦ « CASE ,E12.5,27H * * ^ STREAMWI 
ISE DIMENcion ,e12*5,13H FEET * * * ^ /) 

■^771 FORMAT! IDSEIC.3) 

qoOl F0 PmaT(I3,7E10.3) 

J=MQp( ITFm, 2)4-1 
IF(MWE)1C1,154,154 
ICl CALL PECASE 
NBT=3 
NRT2=Q 
IS=1 
1-^=1 

IF(KP( D) 104,104,103 

103 CALL LINMAT 

104 KR17=KR( 17) 

154 IF( IS+ITFM-2) 153,153,1572 
1572 IF(NITEM-1) 1574, 1262 , 1574 
1574 WRITE (NBT) HI STl .HI ST2 ,HI ST3 ,H I ST4 ,VMAT 
IFf item-1) 157,156,157 

156 IDUM=N6T 
NBT=NBT2 
N8T2=IDUM 
REWIND NBT 
REWIND NBT2 
GO TO 155 

157 IF( IS-1 ) 153.153,155 

1^5 READ! NBT2)HI STl ,HI ST2 ,HI ST3 ,HIST4,VMAT 
GO TO 1262 
153 CONTINUE 

IF! IT-1) 105,105,106 

INITIAL GUESSES FQP PRINCIPAL DEPENDENT VARIABLES. 
CALCULATE (KR(2)=0), INPUT (KR(2)=1), OR USE VALUES 
QPOM PREVIOUS CASE (KR(2)=2). NOTE: LATTER REQUIRES 
SAME eta value. ITS UTILITY IS FOR REPEATED SIMILARITY 
SOLUTIONS. IT OBVIOUSLY CANNOT BE USED FOR FIRST CASE. 

105 CALL FIRQTG 
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ire. CALL R'^FCCN 
KP( 2)=2 
IS=1 

1262 IP(TIMF( n )1053, 1054, 105^ 

11^3 TIMP( l)=-TIMP(n 

WPTTF(KGUT,5)TTME( ITFV) ,S(TS) 

TIMC( 1 )=-t!MF(1 ) 

GO TO l?f 

10*^4 WP!TF(K0UT,3)TIMP( T^rv) ,s( IF) 
WPTTE(^,,3)TT^^P( IT^M) ,S(IS) 
r COMPUTE HTCTGPIC INPPPMftTIQN 

126 CALL HIFTXI 

MWF = C 

C START OF I’^FRATICN LOCP 

1‘^a ITC^rO 

KR< 17)=KP1'^ 

1*^5 RF-rUPN 
FNin 



SUBOOUTTN'E RFC4^F 
PT^AL *^0 

common /FHGC CM /PFM('^C) ,P'^F(53) ,DIJFS ,0 5(50 ,5'.} ), 05(50,^0) 
l,PMn,Mtl,T,0<^IP( 50) ,IPSTP 

COMMON /ETAC CM/E Ti ( 25 ) ,nET^(24) ,050(24) , OflJ ( 24 ) , B 1 ( 2^ ), 
1R2( 24) ,L/P( 100) ,BAir^3,25) 

CnMMON/INTfr'v/KP (2'') , KT N , K GUT , m i ,matij,NFT A, I, I$,N'S 
1,T'^,^^^TME ,N^ M,\»l fo,NML^O,NPNL,!TS,C A~E( 1*^ ), B(S) ,MWE,NO 
2N,ITFM,MTEm,kp 17 ,NeT,NBT2 , 1 GENT ,U MFS , PRr 0, EPS 
COMMON /PPMCPM/TI ME (50) ,*5RE (50 ) , PTET,S (50 ) , POKAo ( 50 ) , PN 
inSE,VK'^P,NOT PC , ini SC (50) ,NSn(lQ ) ,MSO( 1C ) , IPPE,PAONO,CO 
2NE,DEL^, CYCLE 
COMMCN/PRPCCM/PP ,SC ,XM(5 ) 

COMMON /WALCnM/FW(^r)) ,P«nvW(50 ) , TEW, IRHOVW, JPHCVW 

1 E0RM,"T(1H1///////15X4ohmnSTEADY BOUNOAdy layer INTEGRA 
IL matrix PPGCEOUPR // 

12‘^X25H0rPA Rtment OP t:o hnAMTI /2“^X 25HN^ V AL POSTGRAOU 
2ATC ctmcpl /30X2CKM0MTPPEY , ''ALTEQOMtA / / 32X 1 6H I NV ES T I 
3GATICN nr //10X15A4, //02X15 

4HC0NTRCL NUMBERS /11X58H1 2 3 4 5 6 5 9 10 1 

51 12 13 14 1^^ 16 17 18 19 20 /Q X20 1 3// 35X lOHET A VALUES 
6 /4( 5X1P7E10.3/) /) 

2 FOpMAT( 20X17HN0SE RAOIUS 

3 FORMA T( 20X17HCQNE HALP ANGLE 
EOPMAK 2CX17HDFMSITY 

«=: =r>PMAT(2CX17HVI^CnST'"Y 
10 ) 

A P0Rm 4T( 2CX1'^HTI GRin 
1 ) 

7 EOP.ma T( 2CX17HC YC LE TI mc 
I ) 

o FORMA T( IMl //////////31 XIPHSTATI CN TNPIJT PAT A //IXIIHCI 
IETAMCE ,ET ,12X5(1X1 tE10.3) /9 (24XP (1X1PE10.3)/)) 

9 F0PMAT( 1X5HPCKA P , 1BX‘== ( 1 XI ^ElO. 3 ) /9 ( 24X5 (1X10E10,3)/)) 

10 POPMft T( 1 X20HPRE5«^URE patio, MEAN ,3X/p ( 24X5 ( IXIOEIO .3 ) 
15( IXIPEIC.3) /) ) 

11 forma T( / 5X68HTEMPFPA TUPc tttaL mt^n 

1 ocriL. CSCTL. /21 X53HPPESSURE VELG 

2CITY amplitude frequency /7X66H(0EG R) (A 

3TM) (PT/SFC) PARAMETER (HERTZ) /5(5X 

41PF10.3) ///) 

12 FnPMAT( 1X21HFNTO0PY DROP,B'^U/LB R , 2 ) /p ( 24X5 ( IX I PE 10 . 3 
1X5( 1X1PF1C.3)/) ) 

13 fopmaT( 1X22HMASS FLUX,LB/SEC ET**2 , IX'^ ( IX 1 PF 10 . 3 ) /<= ( 2 
14X5( 1X1PF10.3) /) ) 

IF (ITEM-1)3C2,302,1P4 
3C2 CALL INPUT 

WP!TE(KrUT,l)CA SE,KP , (FTA( JN) ,JN=1 ,NETA) 

wpite(kcut,ii)t,ptft,umfs,eps,preo 

NTTEM=NTIMc 

MTIMF=1 

COMPUTE INFORMATION NFcded TQ CCNSIOER DISCONTINUITIES 



=10.5,13H FT 
f=10.5 ,13H DEG 

,13H lb/ft CU 
= lf^E10.3,13H LB SEC/PT 

= 1'5*=10.3,13H EEC 

= IPF10.3,13H 
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303 J=1 
102 1)=1 
inisc( i)=o 

?{ 1)=ABSC ?;(!)) 

IF (NS-1) 105.105,1021 
1021 on 101 1 ^= 2 , NS 
ini5C( IF)=0 
IF(S(IS)) 103,101,101 
103 NSD( J)=IS-MSD( J)+l 
S( !S)=-S( l^) 

IDTSC( IS)=1 
J=J+1 

MSn( j)=is 
101 CONTINUE 
105 NSD( J)=NS-M50( J)+l 
ND!SC=J-1 

IP( !A0S(KR( 6) -2) -1 >207,208,207 
207 IF(NS-1)234, 234, 20*^1 
20-^1 IIS = 2 

LN7^l 

TF(RCKAP(1>) 223,226,226 

223 RAnNO=-POKAP( 1 > 

ROKAPl 1)=0, 

on 229 IS=2,NS 
IF(PnK/5P( IS) ) 224,224,225 

224 TF(KR(6)> 221.221,222 

222 ROKAP( IS)=S( IS) *SI N ( p ADNO/57 . 29 573 ) 

GO TO 22Q 

221 ROKAP( TS)=RA0Nn«SIN(S(IS)/9A0N0) 

229 continue 
GO TC 234 

225 !F (IS-NS) 2251,234,234 
2251 IIS=IS+1 

LN7=I S 

226 00233 !F=IIS,NP 
IF(R0K&P( IS) f 233,233,227 

227 IF(IS-1-LN7) 232,232,228 

228 LN7=LN7+1 

ROKAP( LNZ)=ROKAP( LN7-1 )♦($ ILNZ) (LN7-1 ) ) / ( S ( IS ) -S ( LN Z 
1-1) l^'IROKAPf IS)-ROKAP(LNZ-l) ) 

GO TO 227 

232 I N7=I S 

233 CONTINUE 

234 VKAP=1. 

GO TO 210 

2C8 no 209 IS=1,NS 
20P ROKAP( IP)=1. 

VKAP=0. 

21C CONTINUE 

181 <^TFF = .481E-12 
IF(KR(6)-1)1P2, 182,183 

182 WRITE(KCUT,2)PA0N0 
WR!TP(KrUT,3)CnNE 

103 WRITE! KCUT,4)RH0 
WRITE! KCUT,5)MU 
WRITE! KrUT,6)DELT 
WRITP! kCUT,7)CYCLE 
WR T TE ! KOUT , 8 ) ! S ! I ) , I = 1 , NS ) 

WRITE!KCUT,9)!R0KAP!I ) ,1=1 ,NS) 

WRITE! KCUT ,1C) ! PPF ! I ) , I =1 . NS ) 

WRITE !KfUT, 12) ! 0^1 P ! I ) , I =1 , NS ) 

WRITE! KCUT, 13) ! RHOVW! I ) ,1=1 , NS) 

184 RETURN 
END 



SUBROUTINE INPUT 
PFAL MU 

COMMON /EDGC OM /PE M! 50 , PTE! 50) ,OUFS ,UE ! PQ ,50 ) , PP!50,50) 
1,RH0,MU,T.DSIP! 50) , I OSIP 

C0MM0N/ETAC0M/ETA!25) ,0ETA!24) , OSO! 24 ) , DCU ! 24 ) , B 1 ! 24 ) , 
IB2!24) , LAP! 100) , BA 1! 73,28) 
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^P^^M^^!/INTCnM/KR (20) ,KTN,KnU'r,MATU ,m^T1J.NFTA, I, is, ns 
1,!T,NTIMF ,NAM,NIF0,NNLFQ,NRNL,I TS,CA*^E(15 ) , B ( 8 ) , MW F , NO 
2N,TTF^',NTTFM,KP17 ,NBT,MBT2 , 1 OEMT ,IJMFS • PRFO, F^S 
common /FRMCQ y/TlMF (50) ,''RE(50) ,PTET,^(F0 ) . R OK AO ( 50 ) , dn 
inSF ,VKAF ,NDT SC ,101 SC (50) ,NSO( 10 ) ,M$D(1C ) , IPRE ,R A ONO, CO 
2NF,0FLT,CYCLR 

COMMON /VAR C CM /F (5,25) ,FTM1 (5,25 ) ,FTM2 (*= ,25 ) 
common /WALCOM/FW( 50) ,RH0VW(*^0) , IFW , I RHOVW , JRMQVW 

1 FORMAT(20I 1 ,15AA) 

2 «=npMAT(2I3,F10.^ ) 

3 FORMAK 8F10.4) 

4 format( 5(2X514.7)) 

RP/n (KIN, DKR, CASE 
PP^D (KIN,2)MTMF 



100 

K1 



102 

1C3 

104 

105 
1C6 
ic-' 

1C“^1 

110 

112 

113 



RPAH (KIN,2)NS 

RPAD ( KIN,3) ( S( I P) ,I S=1 ,NS ) 

TP (KP ( 1 )) 103,103,100 

IF(KP( 12))101, 101,102 

RFAH (KIN,?)NETA 

PF/0 ( KIN,3) (ETAd ) ,1=1 ,NF'ta) 

r,n -p 1C3 

CALL 5T0ETA 

IP ( »^R( 6)-l) 10^,104,10^ 

PEAO ( KTN,3)rONF ,PNOPE 

I F( TABS ( KR( 6)~2 )“1 ) 106,107 ,106 

REAH (KIN,3)(RCKAP(I'^) ,IS = l,NS) 

REAP (KIN,3)PTFT 

I«= (KO(2)-l)112,lin,112 

RE/n (KIN,4) ((F(I ,J) ,J = 1 ,NETA) ,1=1 

IF( ITCv-i) 115,113,11^ 

insip=i 



) 



TPPF=1 

no 114 1 = 1 , NS 

114 DP!P( I )=C. 

115 !F(KR(5)-^)118,116,118 

116 IF( TnSTP-ITPV)liP,n7,ii« 

117 PF^n (ki\,3)(D^IP(I) ,I=1,MS) 
PC»r) (KTN,2)^^'^IP 

IIB IF ( IppF-I^fm)12C,11®,120 

119 RFAD ( KIM,4) (PPF (I ) ,T =1 ,NS) 
PFAD (KIN,2)IPPF 

120 JPHOVW=C 

IF( ITCm-1) 121,121,122 

121 IPW=1 



T0H0VW=1 

122 READ (KIN, 3)-^ 

12<= TFf ITFM-1 ) 130,130,131 

130 IF(KP(P)-1)131,135,138 

131 IF( IFW-!-»'EM) 132,133,132 

132 IF( TFW-1) 134,13*^.134 

133 REAP (KIN, 3) (F W( I ) ,I =1 ,NS) 
REAP (KTN,2)IFW 

134 GO TO 13« 

135 IF( IRHOVW-ITEM) 136 ,137,136 

136 GO Tn 139 

137 READ (KIN, 3) (PHOVWn ) ,1=1 ,NS) 
RFAO ( KIN,2)IRH0VW 

138 RFAO( KIN,3)UMF*: ,FPEO,*='PS 
JRH0VW=C 

IF (FRF0)141, 141,139 
13P CYCLF=l./FPco 

BAG=FLOAT( NTIMF) 

oflt=c yclf/bag 

TIMF( 1)=0. 
no 140 I=2,NTIMf 

140 TTMF( I)=TIME(I-l)*^nELT 
GO TO 142 

141 TIMF( 1)=C. 

OELT=0. 

CYCLE=C. 

142 continue 

IF(KR( 16) ) 143,145,143 
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1A3 DO 144 1=2, NS 

144 PPP( I)=PPE(I )*PH0*UMFS^I|MFS/136175.5 7>1. 

145 RETURN 
END 



SUPRnUTINE STDETA 

COMMON/FTACHm/et/^ (25) tOETA(24) , OSO ( 24 ) » DCU ( 24 ) , 
1B2( 24) ,LAR(1C0) ,BAl(T3t2P) 

COMMON/INTCnM/KP(2 0) » KI N tKGUT » M ATI I , MAT IJ , NET A, 
It IT, NT IMP ,NA M,NLEOtNNLEOtNRNL,ITS,CASE{ 15 ), B(al 
2N, ITEM,NITEM,KR17tNBT tNRT2 , I DENT ,U MFS , EREO, EPS 
IF(KR( 12)-2)4,5,1 

1 IF(KR( 12)-4)6,7,2 

2 IF( KP ( 12)-6) 8,9,3 

3 IF(KR( 12)-8) 10,11,12 

4 NFTA = 8 
ETA( 1)=C.C 
ETA(2)=C.5 
FTA( 3)=1.C 
FTA(4)=1. 5 
ETA( 5)=2*0 
ETA(6)=3.C 
ETA( 7)=5.0 
pta( 8)=6.C 
GG TO 13 

5 NETA=8 

ETA( 1)=0.C0 
FTA( 2)=C.90 
FTA( 3)=1.60 
ETA( 4)=2.40 
FTA( 5)=3.20 
ETA( 6)=4.00 
ETA( 7)=5.00 
FTA( 8)=6.0 
GO TO 13 

6 NETA=10 
FTA( 1)=C.00 
pta( 2)=C.60 
ETA( 3)=1.2C 
FTA( 4)=1.80 
Em 5) = 2.40 
ETA( 6)=3.C0 
ETA( 7)=3.60 
PTA( 8)=4.20 
ETA( 9)=5.C0 
FTA( 10)=6.0C 
GO TO 13 

T NFta=12 

ETA( 1)=0.00 
ETA( 2)=C.50 
FTA( 3)=1.00 
FTA( 4)=1.50 
ETA( 5)=2.00 
F^A( 6)=2.5C 
ETA( 7)=3.00 
FTA( 8)=3.50 
ETA( 9)=4.00 
FTA( 10)=4. 50 
FTA( 11)=5.00 
= TA( 12)=6.C0 
GO TO 13 

8 MFTA=14 

FTA( 1)=0.00 
ETA( 2)=0. 40 
ETA( 3)=0.80 
FTA( 4)=1.20 
ETA( 5)=1.60 
FTA( 5)=2.00 
FTA( 7)=2.40 
ETA( 8)=2*80 
ETA( 9)=3.20 



BK24), 

I, IS, NS 
,MWE,NO 
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FT/^( 10) = 3.60 
f-TAC 11)=A.00 
FtA( 12)=^. 50 
O'TAJ 13) = 5.C0 
FTA( 1A)=6.0C 
on in 13 
c mcta = i*2 

FTA ( 1 )=C.OO 
FTA ( 2)=C.2C 
FTA( 3)^0. 40 

( M=C. 60 

?ta( 5 )=r,ec 

( A) = 1.C0 
PTA(?)=1.20 
FTA( ^3) = 1.4C 
6TA( «)=1.6C 
F^A{ 10)=2.0C 
=TAf 11)=2. 50 
FTa ( 12)^3.25 
FT/ ( 13)=^. CO 
PT/ ( 14 )=F .00 
>^T/ ( 1F) = 6.C0 
GO TO 13 



1C NFTA 


= 1P 


FTA ( 


n=c.oo 


FTA ( 


2)=C.3C 


FT/ ( 


3)=C.6C 


F-A ( 


4)=O.FC 


FT/ ( 


*^) = 1.20 


FTA( 


^)=1. 50 


FTA( 


7)=1. 50 


FT/ ( 


8) = 2. 10 


F ta ( 


9)=2.40 


FT/ ( 10 ) = 2. 70 



11 



12 



FT/( 11)=3.00 
ET/f 12)=3.30 
F^Af 13)=3.6C 
FT/ f 14)=3. FC 
FT/f 1*^)=4.20 
cjA( 16)=4.60 
ETA( 17)=*=^. CO 
^TA ( 1 R)=A .00 
on TO 13 • 
NFTa=22 
FTA( 1)=C.000 
2)=C. 125 
FT/ ( 3)=C,250 
^ta(^-)=C. 3 T- 
F^^A ( 5 )=C, 500 
FTA(6 )=C.625 
pta( 7 )=F .760 
FTa ( o ) =C. 

FTA( F)=i.0C0 
ETA ( 10 )=1. 125 
ETA( 11) = 1.25C 
FTA( 12) = 1 .3T*^ 
FT/,( 13) = 1. 5CC 
FTA( 14) = 1.-’5C 
ET^( 15)=2.00C 
FTA ( 16)=2.25C 
p-A i 17)=2. 5CC 
FTAf 10)=3.COC 
FTA ( io)=:3, Foc 

Fta( 20 )= 4 .CCC 
ETA( 2l)=^.0CC 
fta( 22)=5.00C 
GO TO 13 
NFTA=25 
FTA( I ) =0.00 
FTA( 2)=C.20 
*^T/( 3)=C.40 

FTA( 4>=C.60 
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5 ) = 

?TA( 6 ) 
ETA( 7) = 
ETA( = 
FTA( 9) = 
ETA( 10) = 
FTA( 11) = 
ETA( 12) = 
ETA( 13) = 
'=TA( 1A) = 
ETA( 1E) = 
ETA( 16) = 
ETA ( 17) = 
ETA( 18) = 
ETA( 19): 
ETA( 20 = 
ETA(21) = 
ETA( 22) = 
ETA(23) = 
ETA( 2A) = 
^ta(25) = 
13 RETURN 
END 



:C.80 
1.00 
= 1.20 
= 1.40 
1.60 
= 1.80 
= 2.00 
2.20 



40 

60 

80 

00 

20 



= 3.40 
= 3.60 
3. 8C 
=4.00 
4.25 
= 4.50 
= 5.00 
= 6. CO 



SUBROUTINE LTNMAT 

COMMON /ETACQm/ETA (25) ,OETA ( 24 ) , DSO ( 24 ) , DCU ( 24 ) , B 1 ( 24 ) , 
1B2(24) ,LAR( 100) ,RA 1(73,28) 

C0MM0N/INTC0m/KR(20) ,kin,k out, math ,matij,neta. I, is, ns 
1, ! T,NTT ME ,NAM,NLE0, NNLEO.NRNLfTTS, CASE (15), B( 8) , MWE,NO 
2N.ITEM,NITEM,KP17 ,NBT,NBT2 ,I0ENT,UMFS,FRE0, EPS 
COMMON /TEMC0M/SPDUM2 (10) ,DER (50 ) ,DUMM1 ( 25 ) ,SLOPE ( 25 ), R 
1E0UM( 25) ,S0UM1 (50) ,SDUM2(50) ,FWOUM(50 ) , X I CON ( 50 ) , FW CO 
2N( 50) ,DUDS(5C) 

DO 104 I=2,NETA 
DETA(I-1)=ETA(I )-ETA(T-l) 

DSO( I-1)=DE'"A( I-l) *DFTA(I-1) 

Bl( I-1)=B( 3)^DS0(I-1) 

B2(I-1)=2.=«'B1(I-1) 

104 0CU(I-1)=^DETA(I-1)*DS0(I-1) 

MAT1I=3*NETA-2 

matu = nftA+2 

no 105 1 = 1, MATH 
DO 105 J=1,MAT1J 
1C5 BAK 1, J) = 0. 

BAK 1,1) = 1. 

BA1(MATH,2)=-1. 

DO 106 I=2,NETA 
IF (KR(IO)) 1041,1041,1042 

1042 (I-NETA) 1044,1041,1044 
1041 DUM1=B(5) 

DUM2=B(4) 

DUM3=B(3) 

DUM4=B(2) 

0UM5=B( 1) 

0UM6=B( 1) 

GO TO 1043 
1044 0UM1=C. 

DUM2=B( 3) 

DUM3=0. 

DUM4=B( 1) 

OUM5=0. 

0UM6=1. 

1043 BAl ( I-1,I+1)=DUM2*0CU(I-1) 

BAl ( I-l,I-»-2)=DUMl*nCU(I-H 
J=H-META-1 

BAl ( J,I+1)=DUM4^DS0(T-1) 

BAl ( J,I-*-2)=DUM3*DS0(T-l) 

J=J+META-1 

BAl ( J,H1)=0UM64DETA(I-1) 

106 BAl ( J,I+2)=DUM5*DETA(I-1) 

9060 FORMAT(2X1P12E10.3) 
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Q 9 f 4 '"(?X 37 HlT^'F^*P matrix PpR MOMENTUM ECU AT I CNS ♦ I 3 , 2 H 

1 XT 3 , 2 '^H, BEFOPE AND AFTER SOHITTGN) 

IP(KP( 15 )- 1 ) 9 C 62 , 9062,«5061 
Q 061 WRTTP( KPUT,QC 6 ^) maTIT ,MAT 1 J 
on 9063 1=1 , MATH 

906 3 WRT TF( KrUT, 9 C 6 C ) (PAl ( I , J) , J = 1 ,MAT 1 J) 

9062 no 107 LI= 1 ,MATij 
107 CALL MATnMcfRAl(l,LT )) 

IF (KR( 1 F)- 1 ) 9072 , 0072, 9066 
0066 m 9067 1 = 1 , MATH 

9067 WRI TF( KPUT, 9 C 6 C ) ( B A 1 (I , J ) , J= 1 , M AT IJ ) 

9072 NLFQ=MAt 1 I 
NMLEO=M/ TIJ 

nam=nnlfo-i 

NRNL =1 

9902 CPNTINUF 
Return 
END 



SUBROUTINE matONF(E) 
niMFNSICN 1) 

CPMMON/ETACOM/ETA (?*=) ,CFTA(2i^) ,DS0(24) , DCU ( 24 ) , R 1 ( 24 ) t 
1B2( 24) ,LAR ( 100) ,BA1(73 ,28) 

COMMPN/INTCPM/KR (20) ,KIN,K CUT, MATH ,MAT1J,NE'^A, I, IS, NS 
1,T*^,NTIMF,N6M,nlF 0,NNL»=0,NRNL ,ITS,CASE(15),B(8) ,MWE,NO 
2N, TTF v,NI T^M,KRH ,NBT,NRT2 , 1 DENT ,IJMF$ , FPEO, F^S 
I=MAT1I 

DO 201 J=3,NETA 

1=1-1 

2C1 E( I )=F( I)+F( I+l ) 

I=NFTA 

M=NPTA+NETA-1 
K = C 

no 202 J=2,NETA 
1 = 1 + 1 
m = m-H 
K = K + 1 

2C2 E( I ) = f=f T-l)-E( I )+pcTA ( K) *F (M) 

20 3 E( 1 )=-F(l)+nETA (1 ) *F ( K + 1 ) + nSQ(l ) /2*«E(I + 1 ) 

1 = 1 + 1 
K=K + 1 
M=i 

no 204 J=3,NETA 

! = ! + ! 

K = K+1 
m = m+1 

204 '^(M) = E(M-l)-H(M)+0FTA(M)*RfK)+nE0(M)/2.«F(!) 

RR'^URN 

CNjn 



C 



SUBRPL’TIME pFPAY{M,C,C,NNN,IS ,ND) 

nr mens TEN nCNn,l) ,snc2® ) ,C (Nn,l) ,L( 2 °) ,S( 2 B),LL( 2 B),LL 
1L( 28) 

C9MMCN/INTC0M/KR(20) ,KJN,KOUT 

nTPEET INVERSION PPPOEDURF — C IS RFPLACEO BY C*«-l 
N1=N+1 



1 = 1 , 

= I 

113 L( n = i 
11 rONTINUF 



DP 11 
LLLfl ) 



IX=-1 



IF( IS + 2) 111 ,109,111 
1C6 FHPMATdlH LCD ,1=1, 13, 5X (3013)) 

107 F^pmah lEH ( (C ( I , J) , J = 1 ,I3,12H) ,(n( J) , J=l, 13, 
1,13, 15H) BEFORE REPAY) 

ICB FORMAT( 2X1P11E10,3/(12X1R10»=10.3) ) 

109 WR ITC( kOUT,1C7) N^NNN.N 

WPITF( KrUT,lC 6 )N,( L( T ) , 1 = 1 ,N) 



I X=0 



no 110 1 = 1 , N 



6H), T=1 



79 



lie WPITF( KOUT,10B) (cn » J) »J=1 »N) , ( D ( I f J ) , J= 1 , NNN ) 
111 TS=-1 

TRIANGULATE MATRIX 
on 15 I=ltN 
DO 16C M=1,N 

160 S(M)=ABS(C( I ,M) ) 

IP(I.S) 18,16tl6 

IR IS=0 

GO TO 12 

REDUCE ROW I BY PRECEEDING ROWS 

16 DO 17 J=2,I 
K=L( J-1) 

DIV=-Cn,K) 

IF(OIV) 161,17,161 

161 C( I,K)=0. 

D0162 M=1.N 
DIVC=DIV*C( J-1,M) 

S(M)=AMAX1( S(M) ,ABS(DIVC) ) 

162 C( I ,M)=C( I ,M)-i-DI VC 
IF(NNN) 17,17,163 

163 m 164 M=1,NNN 

16A D( !,M)=0( T,m)+0IV*D( J-1,M) 

17 CDNTINUF 

SEEK MAXIMUM PIVOT 

12 DIV=0. 

DO 13 JJ=I,N 
M=L( JJ) 

IF(ABS (C( I,M) )«0IV)13,13,12l 
121 DIV=ABS CC( I ,M) ) 

K = M 
J=JJ 

13 CONTINUE 

SD( I )=DIV/S( K) 

U J)=L( I ) 

L( I )=K 

IF( SD( I )-l.F-8) 130,130,14 

130 C( I,K)=0. 

IF(SD(in 131,131,12 
SINGULAR MATRIX RETURN 

131 IS=-1 

WRITE( KCUT,151) (II ,L(II ) ,SD(II ) ,I 1=1, I ) 

RETURN 

14 DIV=C( I,K) 

C( I,K)=1.0 
K=LLL( J) 

LLL( J)=LLL(I ) 

LLL( n=K 

LL( K)=I 
NOPMALIZE ROW 
IF(NNN) 143,143,141 

141 DO 142 J=1,NNN 

142 0( I , J)=D( I , J) /DI V 

143 DO 15 J=1,N 

15 C ( 1 1 J ) =C ( I , J ) /O I V 
IF(fx) 152,150,152 

151 F0PMAT(24H oiVOT ROW/COL/RES. RATIO 5(I4,1H/I3, 
IIH, )) 

IFC WRITE(KCUT,151) ( I , L (I ) , SO ( I) , I = 1 , N ) 
DIAGONALIZE MATRIX 

152 NM=N-1 

DO 20 1=1, NM 
K=L( I-H) 

DO 20 J=1,I 
DIV=-C( J,K) 

IF(OIV) 19,20,19 

19 C(J,K)=C. 

IF(NNN) 191,191,192 

192 DO 193 M=1,NNN 

193 D( J,M)=D( J,M)+DI V*D(I-H ,M) 

191 DO 201 M=1,N 

201 C( J,M)=C( J,M)^DIV*C(I+1,M) 

20 CONTINUE 



1H/E9.2, 
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c 



r 



TMTfPCH/NGP COLUMNS 
on 30 TI=1,M 



T = ir 

21 J = L( I ) 

L( I ) = T 

IF( J-T )22, 30,22 

22 !F( I«^)2f ,23,25 

23 nn 2 ^ 

2A C(M,I )=C(M,J) 



I=J 

r,n TO 21 

25 IPf I^-J)26,?P,2ft 
2(^ nn 27 M=1,N 
27 rtM,T >=C( w,J) 

T=J 



GO Tn 21 



2P nn 2^ w=l,N 

29 r(M,I )=S(M) 

TS = 0 

30 rnNT|Njyp 
IN'^rPCM^ N'OF DflWS 
nn 40 11=1, N 

I = TI 

31 J = LU I ) 

LL( T )=I 

J-T)32,40,32 

32 IF( IC)35,33,35 

33 no 34 M=i,N 
S ( ) =C f I , M ) 

34 c( I ,M)=r { j,M) 

TFfNNN) 343,343,341 

341 on 34? M=l,^)^N 
^n( M )=p( T ,M) 

342 D( ! ,M)=0( J,M) 

343 IS=I 

!=J 

r,n TP 31 

35 IP( IS-J)36,36,36 

35 nn 37 M=i,N 

37 Cn,M)=C(j,M) 

irnsjMN) 3"^3,373,371 
3“’l nn 372 y=l,NNN 

3“^2 n( T ,M)=P{ J,M) 

373 T=J 

no TO 31 
pn 3Q V = 1 

39 C( !,»^)=S(M) 

IP(NNN) 3°3,393,391 

3^1 no 392 

3C2 n( T ^M)^cp( M) 

3<=3 IS = 0 

40 CONTINUF 

IF(IX) 411,409,411 

407 ^ORM/iT(15H ((C(I,J) ,J = 1,I3,12H),(0(J),J = 1,I3, 
1,13, 15H) 4FTFP RFPAY ) 

409 WRITF( KGUT,4C7) N,NNN,N 
no 410 1=1, N 

410 WR!TF( KOnT, 1C8) (Cn , J) , J = 1 ,N) , ( H ( I , J ) , J = 1 , NNN ) 

411 PP'^UPN 
FND 



6H), 1=1 



SUBPOUTTNF FIRSTG 

COMMON /FTACOM/FTA( 25) ,DFTA f 24 ) , DSO ( 24 ) , OCU ( 24 ) , B 1 ( 24 ) , 
1B2( 24) ,LAR(100) ,BA1(73,28) 

COMMON/INTCOM/KP (2 0 ) , KI N , KQUT , M ATI I , MAT IJ , NET A, I, IS, NS 
1 ,!T,MTIME ,NAV,NLE 0 ,NNLE 0 ,NRNL,ITS,CASE( 15 ),Bf 8 ),MWE,NO 
2 N,ITFM ,NITEM,KR 17 ,NBT,N 0 T 2 , T DENT , 1 ) MFS , F RFO, EPS 
CnMMON/PPMCOM/TIMF ( 50 ) , pRE ( 50 ) , PTE^ , S ( 50 ) , ROK AP ( 50 ) , RN 
lOSF ,VKAP,NDT *^C , I DI SC ( 50 ) , NSD UC ) , MS D (ID ) , I pRF , R A DNO, CO 
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2NF,DELT, CYCLE 

C0MM0N/VARCGM/F(5,25) .FTMl ( 5 , 25 ) ,ET M2 ( 5 1 25 ) 
IF(KR( 2)-l) 112,152,152 
112 nn 128 I=1,NETA 
f( 1,1 )=0. 

0UM1=ETA( I )-A. 

IFf OU^l) 115,116,116 

115 DUM1=-C.1*0UK1 
OUM2=0.25*ETA(I ) 

F( 2,n=OUM2 

3,1 )=nUMl 
GO TO 114 

116 F(2,I )=1.0 
Ff 3,1 )=C. 

114 F(4,I )=0. 

128 continue 

IF(KR? 20) ) 152,152,200 

200 DO 208 I=1,NETA 
IF(ETA( I )-5* )201,201 ,202 

201 F( 1,1 )=ETA(I ) 

F(2,T)=ETA(I)/5, 

Ff 3,I)=*1^(5,-ETA(I) ) 

GO TO 203 

202 F( 1,I)=5. 

F( 2,1 ) = 1. 

F(3,I )=0. 

203 TF(ETA( I )-2* )204,2C4,205 
2C4 F(4,I )=-.l*ETA(I ) 

GO TO 208 

205 IF(FTA(n-5.)206,206,207 

206 F(4,I)=-.067^(5,-ETA(I)) 

GO TO 208 

207 P(4,I )=c. 

208 CONTINUE 
Ff 1,1)=0. 

F( 2, 1)=C. 

Ff 2,META)=1. 

Ff 3,NE'«'A )=C. 

F(4,nF*^A)=0. 

152 CONTINUE 

9<301 F0RMAT(4f 2X7E10.3/) ) 

QQ04 FORMATf 2X19HDEBUG FIRST GUESSES) 

IF(KR(15)) 9902,9903,9902 

9902 rONTTNUF 
WRITEf K0UT,99O4) 

WRITE (KCUT,9901) f (F(I ,J) ,J=1 ,NETA) ,T=1,4) 

9903 CnN'TINUF 
RFTurm 
FND 



SUBROUTINE PcFCON 

PEAL MU 

COMMON /EDGCOM/pfm( 50) , PTE (50) , DUES ,UE ( 50 , 50 ), PE ( 50, 50 ) 
l,PH0,Mll,T.DSTP(50) ,IDSIP 

COMMON/ETACOM/ETA(25) , DETA ( 24 ) , OSO ( 24 ) , DCU ( 24 ) , 0 1 ( 24 ) , 
1B2( 24) , LAP (ICO) ,BA 1(73,28) 

C0MM0N/HISC0M/TZ,T1 ,T2 ,TF(25) ,C 1 ,C2 , C3 • C4 , BET A, ZM(5, 25 
1) ,XI( 50) ,HF( 25,5) ,HUE ,HHUE,DLX2 ,C3M(50 ) , BET AM (50) 

common/intcom/kr(20) ,kin,kout,mati i,matij,nfta, I, IS, ns 

1, IT,NTIMF,NAM,NLEQ,NNLEO,NRNL ,ITS,CASE( 15), 8(8) ,MWE,NO 
2N,ITEM,NITEM,KR17jNBT,NBT2 ,IOENT,UMFS,FREQ, EPS 

C0MM0N/PRMC0M/TIME(50) ,PRE(50) , PTET , S ( 50 ) , ROK AP ( 50 ) , RN 
lOSE ,VKAP,NOISC ,TDI SC (50) ,NS0(10 ) ,MSD(10) , I PRE , RA ONO, CO 
2NE,0ELT,C YCLE 

COMMON /TEMC0M/SP0UM2 (10) .0ER(50).0UMM1 (25 )iSL0PE(25),R 
lEOUMf 25) ,SDUM1 (50) ,SDUM2(50) ,EW0UM(50 ) , X ICON (50 ) , EWCO 
2N( 50) ,DUDS(50) 

COMMON /WALCOM/EW( 50) ,RH0VW(50 ) , lEW , I RHOVW , JRHOVW 

1 FORMAT( 1H1//////////29X21HSTATI0N COMPUTED DATA /30X5H 
lTIMF=r 1PE10.3.4H SEC /) 

2 FOPMAT( /5X68HSTREAMWI SE XI BETA 
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1 FDGF STATIC /5X69H0I MENS I ON (LB/ 

2 Velocity pressure /7x6 

35H(FFET) SEC)*«2 (FT/SFC) 

A (ATM) /) 

3 phrmaT( 5( 5X1FF10.3) ) 
ooQl '^nRM/^T<E(5X10E10.3)) 

1=:32 IFf TTEM-i) 605tl538,605 
1S38 IOSIP=l 
I^P,F = 1 

DO l*^AO 1 = 1, NS 
ic;40 D<^!P( I ) = C. 

1=^3^ (KR( 6)-3) 1535,1541 ,1535 
1535 IP (PPF( 1) )154ia542,1541 
1*=A2 DO 15A3 1=2, NS 

TP(PPF(I)) 1544,1543,1544 
15^4 L=! 

GO TO 1545 
15^.3 CONTINUE 

15^5 PMnSP = Sf L) /SCRT(1.-PRE(D) 
no l«^46 1=2, L 

154 6 Doen-i) = i,-s( T-l) /RNCSE^S(I-l) /RNOSE 
1541 no 1531 T=1,NS 
DTF( I )=PTET 

1531 n*=M( T )=PTET«PRF(I) 

8002 '>l=3.14ic<? 

no 8 1 = 1 , NS 

on TI = 1,NITEM 

nDQC=( ptet-PEM( T ) ) /PHC 
0^05 = 136 17‘^.57*DPDS 
UE( I ,I I ) = <^ORT(DPDS) 

T ,I T )=UE( I ,I I )*a . + FPS*COS ( ( FPEO«S( I )/UE( T , I I ) )> 
IFLOATf II) «2. ^PI /FLOAT! NITEM) ) ) 

'5'=( I ,I n=PTFT-( RHP+UE ( I ,I I )** 2 . /1361T5.57 ) 

8 rnMTjMuc: 

2C^ IF(KP(15)) 9503,9904^9903 
9903 CON^INUF 

WRTTF(K0UTta901)UF(I ‘^,ITFM) ,°TP (1$) ,t^rhC,MU 
0904 ''HNTINU'^ 

r FHOC*^? PUDS TO ZFRC PCR FLAT PLATF TEST, 

I*" (KR(7))81,83,«1 
PI nn P2 JS=2,NS 
8 2 UF( J5, I'^EMI^UF! JS-1 ,IT9Mj 
°3 CONTINUF 

C END rip EDGE PROPERTY LOOP, START OF PEta AND XI 

C CALCULATION. 

6C5 XI ( 1)=C. 
f = Nni cr + 1 
P'O 111 II=1,J 
K=Nrp( ID 

f T T ) 

1 * ‘ 

LI. =K + M- 1 

IF (IT-1) 6052, 6C52, 402 
60*^2 IF (KP(6)-1) 40C,401,4C2 
C AXIPYMFTPTC BLUNT 

4CC no 403 I=M,LL 

SnUM2( T)=S(I )^S(I) 

IP( S( n ) 403,403,4031 
4031 PUM1=RCKAP(I )/S(I) 

X!CGN( I ) =UE( I ,I TFM) / S( I ) ♦PHG/4. *mij«oum 1 ♦ PU M 1 
4C3 SPUMK ! ) = FnijY2(I )«SnUM2(I) 
ni=4. 

BET&M( 1) = 0.5 
GO TO AC6 

r o\ ANAR PLUN'^ 

4C1 nn 404 I=v,LL 
SnUM2( T )=S( I ) 

TF(P(I)) 404,404,4041 
4041 XIC0N( I )=UE( I *ITEM)/S(I )*RHG/2.^MU 
404 SniJMK T )=5( I )*S( I) 

01 = ?. 

BFTam( 1 )=1. 
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c 



AG 2 
AG7 



40<5 



Ace 

AC5 



40*^2 

A051 

AC6 



PA 

«*5 

66 

A066 

A061 

r 

A063 



A06A 

^065 

A062 

AlCl 

AIG 



111 

9603 



137 

139 

16P 

139 

”^043 

C 

70^5 



■^049 

2291 

7048 



22® 



C 

c 

232 



GO TO A06 

19 (KR(6)-2) 408,407,408 
AXISYMETPIC SHARP 
DO 40® I=M,UL 
$DUM2(I)=Sn)*S(I) 
r*UMl=ROKAP( I )/S(I) 

XirON( I )=PHO*UF ( I , TTFM) ♦MU^OUMl /3**DUM1 

SDUMK I) = S(I )^SDUM2(I) 

xi( 1) =xicrN(i)^sa)«S(i )♦$(!) 

0T=3. 

IF (II-l) 4051,^051,406 
PLANAR SHARP 
no 405 T=M,LL 
S0UM2( I )=S(I ) 

XICQN( I )=PHO*UF( I ,!TEy)«I^U*ROKAP(I )*ROKAP(I ) 
SDUMK I ) = S(I ) 

01 = 1 . 

IF (II-l) 4052,4052,406 
XT( 1) =XICON(l)*S(l) 

MM = M 



FW( 1)=0. 

CONTINUF 

CALL SLOPQ(K,S( M) , UE ( M , I TF m) , DUDS ( M ) ,OER ( M ) ) 

FORCES DUDS TO ZERO FOR FLA^ PLATE TEST. 

IF (KR(7))64,86,84 
on 85 JS=1,NS 
OUDS( JS)=0.0 
CONTINUE 

IF (KR(6)-2) 4066,4062,4062 
I® ( M-l) 4061,4061 ,4062 
IF (PNOSF) 4063,4063,4064 
DUES COMoyTED BY SLOPO 
DUES=OUPS( 1) 

GO T0 4065 

OUFS FopM effective NCSE RADIU9 USING NEWTONIAN FLOW 
OUFS=SORT( 2. /RHO^PF (1 , I TF M) *32 . 1^405 ^2 1 16 . 2 ) / RNOS E 
XICON( l)=RH0*MU/(2.*VKAP+2 . ) ’^‘DUFS 
CALL SL0P0(K,SDUY1 ( M) ,XICONCM) , DER ( M) , X I ( M) ) 

IF (LL-MM) 111,4101,4101 

no 410 l=mm,ll 

RFtam( L)=2./0I«XI ( L) /UF(L, ITEM) ’^S ( L ) /S DUMl ( L ) ♦DUOS (L )/ 
IXICONI L) 

CONTINUE 

FOPMAT(8E10.4) 

CALCULATION CF C3 MATRIX 
no 138 1=1, NS 
IF (KR(6)-1) 137,137,158 
IF ( I-l) 139,139,156 

C3M(I )=-$ORT(RFTAM(I ) / ( DUE S+RHO ♦MU ) ) 

GO TO 138 

C3M( I )=-S0RT(2.*XT (I) ) /(RHO^ROKAPI I) ♦UE ( I , IT EM )«MU ) 

CONTINtlF 

JOHOVW=1 

IF ( JR HCVW) 7047,7047,7045 

CALCULATE FW IF PHOVW GIVEN 

J=NDISC+1 

00 7046 11=1, J 

K=NSD( II ) 

M=msh( II) 

ll=k+m-i 
00 209 I=M,LL 
IF (KR(B)) 7049,7049,2291 
RHOVWf I )=RHGVW( I ) ♦C3M(I ) 

IF (II-l) 7048,7049, 230 
IF (KR(6)) 229,229,230 

VALID AT AXISYMFTPfC STAGNATION POINT ONLY 
FWCON( T )=-RHCVW(I) /(2.^C3M(I) ) 

IF( I-l) 20®, 209, 232 

MODIFICATION FOR AXISYM£TRIC BLUNT AWAY FROM 
STAGNATION PCINT. 

FWCON( I)=FWCON( I)/S( I)*ROKAP( I) 

GO TO 209 
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VALID FCP ALL PLANAR 
?.30 FWCONK I )=-RHOVW( T ) /C 3 M(I) 

IP(KR( f )- 2 ) 209 , 236 , 20<5 
MnniFTCATION FOP AXI$YMETR!C SHARP 

236 ^WCON( n=cviCON( I)/S( U^RGKAPI I)/ 2 . 

20<2 rONTINUF 

RWnijM( 1 )=FWCCN( 1 )#S( 1 ) 

IF(KP( 6 )- 2 ) 241 , 237 , 2^1 
: MODI^TCATinN FOP AXISYMFTRIC SHARP 

237 cwnUMf 1 )=FWDUM{ 1 )«S( 1 ) 

2^1 continue 

7 C ^6 CALL SLCPO(K,SDUM 2 (M) ,rwCON(M) ,OFR(M),FWDUM(M)) 
on 126 T = 1 ,NS 
IF( T - 1 ) 124 , 12^,123 
124 ic (KR(f)-l) 133 , 133,123 
133 IF ( I ) ) 113 , 113,123 
113 FW( 1 )=RHCVH( 1 ) 

GO ^0 126 

123 FW( I )=PWDUV( I ) /SOPT( 2 . *XI ( I ) ) 

126 continue 

7047 WP!TF(KrUT,l)TTMF(TTFM) 

WPI^FI K 0 UT, 2 ) 
on -^050 != 1 ,NS 

70 ?C WR!TP( KnjT ,3 )S( I ) , XI ( I ) , RE TA M ( I ) ,l)E ( I , IT FM ) , P E ( I , IT EM ) 
FNP 



3 



1 

2 



3 

4 
c; 

6 



subroutine SL0FQ(N,X,Y,S,7) 
PIMEMCIOn X( 1) ,Y(1) ,S(1) ,7(1) 
!F(N-l) 

S( 2)=( Y( 2)-Y(1) ) /( X(2)-X(l) ) 
S( 1) = S(2) 
nc=S( 2) 

no 7 1=1, N 
TP( I-H-N) 2,1 ,6 
OR = OC 

IR (1-2)7, 6, 5 
XOT=x( T )-X( I+l ) 

XTT=x( H-l)-X(! + 2) 

XT(7 = x( !+2)-X(I ) 

AA=Y( I )/( XCT^XTC) 

XnTT=X(?T*XTT 
AB=Y( T+1 ) /XOTT 
AC=Y< !+?)/( XTT^XTP) 

AAA=A A^XTT 
ARP=A3*XT0 



ACC=AC*XC'!- 
OA=OC 
QP=S( T ) 

Or=S( 1+1) 

S( I )=AA«( XT0-X(?T) + ABB-ACC 
S( T+1)=AB«( XCT-XTT)+ACC-AAA 
S( T + 2 )=AC«( XT't-XTO)+AAA-A 0R 
IF( 1-2)7, 5, 4 
<^( T ) = ( I )+QA) /2. 

S( I )=( S( I )+OB) /2. 

XD=X( I )-X( I-l) 

YS=Y( I )+Y( T-1) 

Sn=S( ! )-S( !-l) 

SS = ^( I ) 

7( I ) = 7. ( I-1) + X0/2.«(YS-XD/6.«SD) 

S( I ) = SS 

CON^TNUE 

return 

END 



SUBROUTINE HISTXI 
REAL MU 

dimension FD(4) 

COMMON /EDGCOM/PFM( 50) ,PTE( 50) ,OUES ,UE( 50 ,50 ) , PE( 50, 50) 
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URHn,MU.T,0SIP(50) •lOSTP 

COMMON /FTACnM/FTA( 25) , DFTA ( 2A- ) , OSO ( 24 ) , OCU ( 24 ) , 81(24) , 
1B2( 24) ,LAR(lOO) fBAl(73,28) 

COMMON/HT SCOM/TZ.Tl ,T2,TF(25) tC 1 ,C2»C3,C4, BET A, ZM(5t 25 
1 ) ,XI( 50) ,HF(25,5) ,HUF,HHUE,OLX2 ,C3M(50 ) » BETAM (5o ) 
CnMM0N/INTC0M/KP(2 0) , KI N,K0UT ,MAT1 1 tMATlJ.NET A. I, IS .NS 
l.ITtNTTMF , NAM tNLEO.NNLFO.NRNL, ITS t CASE ( 15 ) , B ( 8 ) . NWE , NO 
2N,ITFM^MITFM,KP17tNBT,NBT2 , 1 DENT ,UMFS . FREQ • E^S 
C0MM0N/PRMC0M/TIME(50) .PRE(50) , PTET , S ( 50 ) . ROK AP ( 50 ) , RN 
lOSFtVKAPtNDI SC, ID! SC (50) ,NSn(10 ) ,MS 0 ( 10 ) , I PPE , RADNO, CO 
2NF,DELTtCYCLE 

C0MM0N/VARC0M/F(5,25) ,FTM1 (5 ,25 ) ,FTM2 ( 5 , 25 ) 

initialize axial variation TERMS 

NUL = 0 

DLX1=0. 

on 120 1=1,4 

120 PD( I )=0. 

IP (KR(6)-1) 400,400,401 

400 IP (TS-3) 154,153,153 

401 IF( IS-?)154,153,153 

1^3 IP(KP(3)-1) 154,143.143 

143 DLX1=AL0G(XI(IS)/XI(IS-1)) 

TP (KR(6)-1) 402,402,403 

402 IP (IS-3) 154,155,15-? 

403 IP( IS-2) 154,155,157 
157 IPfKP(3)-l) 154,155,144 

154 D7=0. 

D1 = C. 

02=0. 

m=npta-i 

no 140 1=1,4 
DO 140 J=1,M 
7M{ I , J)=0.0 
140 continue 

DO 141 I=1,NETA 
00 1^2 J=l,3 
HP( I ,J)=C. 

142 continue 
HF( I ,4 )=0. 

1^1 HP( I,5)=C. 

OLX2=0. 

GO TO 130 

COMPUTE TWC- OP. THREF-PQINT DIFFPRENCE RELATIONS 

155 07=2. /OLXl 
Dl=-D7 
02=0. 

GO TO 145 

144 J=!DISC( IS-1) 

IF( J)121,121,155 

121 DZ=DLXH-DLX2 

01=-DZ /(nLXl*0LX2) *2. 

02=nLXl/(DZ«0LX2)*2. 

0Z=-D1-P2 

145 0LX2=0LX1 

122 FD(3)=Dl«P(4,l)-»-D2=«'HP(l,4) 

DO 147 T=2,NETA 

Pn( 1)=D1«F(2,I )+D2*HF(I ,2) 

FD( 2)=01*P( 3,1 )+D2*HF(I ,3) 

P0( 4)=fD( 3) 

F0( 3)=ni*F(4,I )-»-D2*HF(I ,4) 

147 CALL TAYL0R(0ETA(I-1) ,FD(2) ,FD(1),ZM(1, I-l) ) 

SAVE HISTORIC VALUES 
162 nn 164 I=1,NETA 
HP(I,4)=F(4,I) 

HF( I,5)=D1*F(1,I )+02*HF(I ,1) 
on 164 J=l,3 
HF( I, J)=F( J,I) 

164 CnN-^INUF 

COMPUTF GROUPINGS WHICH DEPEND ON DZ 
13C Cl=l.+D7 

BETA=BE*^AM(IS) 

C2=-2.^(C1+BETA+DZ) 
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C3=r.3V( IC) 

C6=RHn + MU’«'C3*C3 

C CnMonir Twn TR THRFR PPINT OTFFFRPNCF TIME OPRIVATIVFS 

yr(*,R( ^)-l)20C,2nO,201 
2C0 !«=( T'!'Fm-2)202 ,2C3,2CA 

201 Tf^UTPM-2)202,2G3,203 

202 *^7 = 0. 

T1 = 0. 

T2=0. 

no TO 2CF 

203 T7=l./prLT 
t1=-T7 
T2=0. 

on TO ?.C^ 

20'^ T7 = 1.‘^/DFLT 

Tl=-2. 

72=.5/PFLT 

C 3AVF HI ATOPIC VAIIJFR, 

2C= Pn 20^ 

I) =T1+Ftmi( 1 ,i ) +T2«FTM2 ( 1 ,1) 

206 CONTTNIT 
20“^ coNTiNjur 

9opz. FORMA T( ^XT2,6XI2/( RX1O10F10.3) ) 

K’P ( 17) ) OFO? ,9006 ♦^^0'=^ 
co^JTT^MJr 
WPTTPf vruT,OF07) 

qoC7 T( 2X2^HPEPUG IF,DLX1 ... 7 v,hf,TF) 

NF*rft Ml=^'^ ta - y 

W0ITF( KfUT,9O0^ ) I S,I TFM,DLX1 ,0LX2t07,01 ,02,T Z,T 1,T2, 

ir. l,C2,C3,r4,Pp, M I , J) , J=1 tNFTAMl ) ,1 = 1,4) , ( (HFf I, j ), 
2J = 1,*^) ,T = l,NrT/i ) 

P^06 CONTINUF 
Rotijraj 
FHP 



OUBPnu'TTMr t^yLGR (0,FV,F,p) 

ph(1) ,p( 1) ,o(i) 

rriMMCM/jNiTCOV/KP (20) , KI N ,K COT , m ati l ,HAT1J,NFT A, I, IS, NS 
1,TT,NTT ,N^N',NLP0,NNLF0,NPNL,I'»’S,CASF(19) ,B(3) ,MWE,MO 

2N,!TFW,Mtcm,kp17,^rt,mbt2 ,IDENT,,UMc<5^FcFq, FPS 

n2=n:«tn 

IF (KP( 10) ) 1,1,2 

2 IP ( I-^>^TA) 

A FP=0. 

P( 1) = ( ( ( PM( 3 ) / A.-FD/2A. ) «D-F(2) 72. ) *D+F( 1 ) )^'0 
P( 2) =( ( ( FP 73C.“P^( 3) /R. ) «n+P(2) 73. ) ( 1 ) 7 2 . )«D2 
Pf 3) = 0. 

0(4) = (( (CV/2C. -FP772. ) *0-F (2) 7R . )’5'0 + F( 1 ) 76 . )*02«D-p (3) 
1(3) 

on TO 3 

1 Pn=R( 3)-cM(3) 

P( 1 ) = ( ( f Pf 3) 76.-FP724. ) *n-F( 2 )/2. ) *0^F( 1 ) )«0 

P(2)=( ( (Pn73C.-F(3)/R. )«P+P(2)73.)’^0-P(1)7 2. )*02 

P( 4)=-( ( (F0/2 52.-F (3) 772.) «P^-F ( 2 ) 730 . ) «D-F ( 1 ) 72A.)=«'P2« 

n 

P(3) = ( ( (P( 3) /20.-FP/72. )’<'D-F(2) 7«. ) *0+P( 1 ) 76 . )*02*D-P ( 
14) 

3 CnM^INUr 
PPTUPN 
EMP 



SURROUTIN'F ITFPAT 

common /PUMP nv/BOMP ,COR WA,EASF, !CGRM,wr)GT , 1777, ISP, IX 
CnMMnN7FRPCDW7PLP(73) ,FLA(73) , FLFM, ELM ( 3 ) , FLMm, IFLM,NE 
UM,TLMM,DPL( 73) ,FNLF (23) ,FNL(2B ) ,pNLEM, pNLM( 3 ) , ENLMM, I 
?pNLM,\'PNLM,i nlmm,PFNL(28) ,HRNL( 6) 

CnMM0N/INTCnM/KP(20) ,KI N,KPUT,MA'»'l I .MATIJ.NFTA, I, I^,NS 
1 , I7 ,NTTMp ,NA v,nlE0,NNLF0,NRNLiT 7’S,CAPF( IP ) ,8( 8) , MWE,NG 
2N,ITPy,NITPM,XRl7 ,N87,PBT2 ,I0ENT,UMFS ,FRPO, FP$ 
CnMMpN/VARC0V/P(5,25) ,FTW1 ( 5 , 2F ) ,FT M2 ( F , ?5 ) 
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lei 



189 

192 

19A 



1Q21 

162 

159 

160 



16C1 



161 

1C3 



FORMAT(22H NCN-C CN VFRGFNT OUTPUT /) 

FORMA T ( 6X1 3, 1 XF 9.3,2 ( 1XF9.4 ) ,2 (2X13, 2X1PF10.3)) 

FORMATf /10X62HTTERATFD VALUED DAMP MAX LI 

INEAR MAX MOMENTUM /6X62HITS TIME FPPW 

2 ERROR ERROR) 

CL0CK3=0. 

ITS=I T5+1 
CALL LINCER 
N0N=0 

CALL NONCFR 
FPPW=F(3,1) 

IF(KR(4)4KR( 17 )+KP( 19) ♦NON) 189, 189, 192 

IF( ITS-1)192,192,194 

WRITE(K0UT,7) 

CL0CK1=C. 

btime=clocki-clockb 

BT!ME=.01*RTTME 

CLOCKB=CLOCK1 

WPITE( K0UT,5)ITS,BTI ME ,F PPW , EAS E , I FL M, E LMM, !FNLM,FNLEM 
IFf ELMM+FNLMM^.. 001 /(EA$F*RUMO)-.00 2) 162, 162, 159 
N0N = 0 
RETURN 

TP(ITS-5C) 161,160,160 
WRTTE( K0UT,1 ) 

WRITE (6,1) 

IF ( ELMM+ENLMM-.C2) 162,162,1601 
N0N=1 

return 

ITERATE OP OUTPUT 

TF(KR(^)) 181,181,193 

M0N=-1 

RETURN 

END 



SUBROUTINE LINCER 

COMMON /ERRCOM/FLF( 73) ,ELA(73) , F LEM, ELM ( 3 ) , ELMM, I ELM, NE 
1LM,ILMM,DFL( 73) ,FMLE (28) ,ENL(2S ) .FNLEM, ENLM( 3 ) , ENLMM, I 
2PNLM.NPNLM.T NLMM,DFNL(28) ,DRNL(6) 

COMMON /ETACOm/FTA( 25) , DETA ( 2A ) , DSQ ( 24 ) , OCU ( 24 ) , B 1 ( 24 ) , 
1B2( 24) ,LAP(1C0) ,BA1 (73 ,28) 

r0MM0N/INTC0M/KR(20) , xf N , K OUT , M j ,MAT1J,NETA, I , IS, NS 
1,IT,NTIMF ,NA M,mLEO,NNLFO,NRNL,IT$,CASE(15),B(8) ,MWF,NO 
2N,TTEM,NITEM,KR17,NBT,NBT2 ,I0ENT,UMFS,FREQ, FPS 
COMMON /PR MCOM/TI me ( 50 ) , PRE ( 50 ) , PTET , S ( 50 ) , ROK AP ( 50 ) , RN 
10SE,VKAp,N 0T5C ,ini SC (50) ,NS0(10) ,MSO( 10) , IPRE,RAONO,CO 
2NF,0ELT,C YCLE 

C0MM0N/TEMC0M/SPDUM2 (10) ,DER(50 ) , OUMMl ( 25 ) .$ LOPE ( 2 5 ) , R 
1FDUM( 25) ,S0UM1 ( 50) , SnUM2 (50) ,FWDUM(50 ) , X I CON ( 50 ) , FWCO 
2N( 50) ,DUnS('^C) 

COMMON/VARCOM/F( 5 ,25) ,FTM1 (5 ,25 ) ,FTM2(5,25) 

: evaluate LINEAR ERRORS FOR MOMENTUM AND ENERGY 

NFLM=2 
NFNLM=NELM 
DO 401 T=2,NETA 
IF(KR(1C)) 4000,^000,4001 
4001 1= (I-NFta) 4003,4000,4003 
^000 DUM1=B(^) 

0UM2=8( 5) 

0UM3=B(2) 

0UM4=B( 3) 

DUM5=B( 1) 

OUM6=1.0 
GO TO 4002 
4003 9UM1=B(3) 

OUM2=0. 

0UM3=B( 1) 

DUM4=0. 

DIJM5=1 .0 
DUM6=C. 

400 2 FLF( T-1)=-(F(1 ,1-1 )+DETA(I-'1)^F(2, I-1)+DS0(I-1)*B(1)*F 
1(3,I-1)>DCU(I-1)*(DUM1#F(4,I-1)^0UM2*F(4,I) )-F(l, I) ) 
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M=I+MFT.6-1 

PLf ( N')=-( P( 2,1-1) ♦DFTA (1-1 ) (3 , 1-l) +0^0( !-!)♦( DUM3*F{ 

14, !-l)+0UM4«P( J,I ) )-F(?,n ) 

M=r! + 2«NFTA-? 

401 FLF( V>)=-(F(3,I-l)-f0ETA( I-l ) ^DUM5*(F(4, I -1 ) +0I I M6* F ( 4 , I ) 
1)-^C3,I)) 

FLF(^!ETA )=-cf2,l) 

C DF'^FFMINF MAXIMUM LINEAR ERRORS 

4C4 CALL ABVAXfM^TlT ,FLF,ELMM,! clH) 

(I form DRCnUC'’’ OF A*«-l ANP LINEAR ERRORS 

46F CALI MATPNF(FLF) 

473 CONTINUE 
RETURN 
ENP 



subroutine ARMAX(^’,X,X^<,I) 
dimension X( 1) 

1 = 1 

XM = ABS (X(D) 

I'^ (N-1) 4,4,S 
E DO 3 J=?,N 
XT=ABR (X(J)) 

TP(XM-XT) 2,3,3 

2 XM^XT 
T = J 

3 CONTINUE 

4 XM=Xf T ) 

return 

FMD 



SUBROUTINE NCNCFR 
PEAL WU 

DIMENSTPN C0RAR(25) 

COMMON/BUMCOM/BUMP ,C0P V4 ,EASE , IC0RM,WD0T , ISO, IX 

CnMMON/FOGCnv/pEMCco) ,ptl(EO) , DUES ,UF { 50 , 50 ) , RE ( 50, *^0) 
l,oMn,My,T,r)sip(:^0) irsiP 

COMM0N/ERRC0M/FLFf73) ,FL4(73) ,F LEM, ELM ( 3 ) , ELMW, IplM,NE 
lLM,TLMM,ncL(73) ,FjMLFf2P) ,FNL(28 ) ,FNLEM, ENLM ( 3 ) , ENLMM, I 
2^NLM,NENLM,T NLMV,DFNL (2 8) ,DRNL(6) 

COMMON /FTACnv/ETA( 25) ,PETA(24) ,DS0(24) , Daj ( 24 ) , B 1 ( 24 ) , 
18 2(24) , LAP (100) ,BA I ( 73,2 8) 

C0MM0N/HTSC0w/T7,ti ,T2 ,tp(25 ) ,C 1 ,C2 , C3 , C-^ , BET A, 7M(5, 25 
1) ,XT( 50 ,HR( 25,5 ) ,HU^ ,HH1JE ,0LX2 ,C3 M ( 50 ) , BET AM ( 5? ) 
common /TNTC CM /KR (2 0 , KT N , KCUT , ^ ATI I , MAT 1J.NET A, I, IS, NS 

I, 1"-, NT I MR ,NAM,MLE0,NMLEQ,NRNL,ITS,CAS6(15 ),B(8) ,MWF,N0 
2N, ITRM,NTTcm,KR17,NBT ,NBT2 , 1 DENT ,1J mf$ , FR FQ, ER$ 

C0MM0N/NCNC0M/A^(2 8,2 8) ,DVNL(2E ) ,C'? , C6 , Ct , CB , C9 , C 10 , C 1 

II, r 12 

COMMON/RRMCGM/’^IME (50) ,PRF(50) , PTRT , S ( 50 ) , ROK AP ( 50 ) , PN 
incp ,VKAF,NOI , IP! SC (50 ,NRO( 10 ) ,MSD( 10 ) , I RR*= , RAPNCl, CO 
2'^'E,0FLT. CYCLE 
COMMPN/PRPCCm/pp,sC , XM(5 ) 

C0MMON/TFMC0M/conuM2 ( 10) ,DER (50 ) ,DUMM1 (25 ) • S LOPF ( 25 ) , R 
lcr>HM( 25) ,SDUM1 ( 50) ,S0UM2( 50) ,FW0UM(50 ) ,XIC0N(5C ),RWCG 
2N( '^O) ,PUPS(5C) 

common /V 4RC0M/F f 5,25 ) ,PTM1 (5 ,25 ) ,RTM2(®= ,25 ) 
CnMMON/U-ALCCM/RW{50) ,RH0VW(50) , !FW, I RHPVW , JPHOVW 
FOUI VALENCE (CORAR(l) ,AM(1) ) 

'^ASF=AMIM(FASP’«'2. ,1. ) 

IF ( T TC - 1) 10,‘=,10 
E^RF^l. 

BUMP = 1.0 
4 IC0RM=1 

COPM4 = 1,= + 10 

10 CONTINUE 
IX = 0 

C EVALUATE coefficients AND ERRORS FOR NONLINEAR 

C EOUATIDNS, 

C INTTIALITP AM MATRIX 

DH 15 I=1,NNLE0 
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no n no 



FNLF( I )=0. 

DO 1*^ J=1,NNLE0 
15 AM{ I , J) = C. 

C EV/L. OPCUPING^: WHICH CHANGE DURING ITERATION BUT ARE 

C NOT F(ETA ). 

IF{ ITEM-2)30,31 ,32 

30 CS=0. 

GO TO 33 

31 C5=TZ*ALCG(UE( I S,ITEM) )^ti ♦AL0G(UE( IS, ITEM-1) ) 

GO TO 33 

32 C5=TZtA!.PG(UE( I S , I TF M) ) +T1 «AL0G1 ) ) +T2^ ALOG (U E ( IS , ITEM- 
1(UF( IS,ITEM-2)) 

33 CONTINUE 

EVAt COEFS, AM, AND ERRORS FOR EDGE BOUN CONDS 

NEXT, evaluate NONLINEAR ERRORS 
45 RNl.EfNFTA4'l)=-(l.-F(2,NETA)) 

*^NLE( NETA<-2)=-F ( 3,NETA) 

PNLE( I )=FW(I S)-F(l ,1) 

NEXT, evaluate original AM COEFFICIENTS 
AM( 1,1)=1. 

AM(NETA+2,2)=-1. 

FNLE(NFTA+1)=FNLE( NFtA+1 )+FLE(2«NETA-1 ) 
no 50 j=i,matij 

CO AM(NFTA+1, J)=AM(NETA-H ,J)+BA1 (2*NETA-1, J) 

START OF major DC LOOP FOR EVAL OF COEFFS AND ERRORS 
AT each STATION. 

DO 120 I=1,NETA 

C TEST TO BYPASS COMMANDS THAT CANNOT BE PERFORMED AT 

C FTACl). 

IF (1-1)75,75,55 
C5 CALL IMONE 

C FVAL GROUPINGS WHICH ARE USED AT l-l AS WELL AS AT I 

75 C6=HF( 1,5) 

C7=F( 3 t )+F{ 2,1 )*(C1*F(1 ,n+C6) 

CP=C4*(tz-C5) 

C9=c iTF( 2 , 1 )-ce 

C10=C6+C1^F( 1,1 ) 

cii=-c4«TFn ) 

C12=C9*FC 1,1 ) 

IF (I-l) 100,105,100 
C RACK TO CONSERVATION EQUATIONS 

ICC CALL lONLY 
1C5 IF (KR(17)) 120,120,115 

110 F0PMAT(21H all THE COEFF IC I ENTS / ( 1 X 1 PI 2E 10. 3 ) ) 

115 WRITF(kOUT,110)C1,C2,C3,C4,C5,C6,C7,C3,C9,C10,C11,C12 
IX=-2 

120 continue 

1^5 IF(KP(19)) 170,190,170 
170 continue 

WR!-^E( K0UT,175) 

175 FORMAK 2X10H0EBUG FNLE ) 

IPO F0Pmat(/(12X1P10F10.3)) 

WR I TF ( KCUT , 1 50 ) ( FNLE ( J ) , J=1 , MAT 1 1) 

IPO continue 

DO 1R<5 III=1,MAT1J 
IPS f!ml( IIT)=FNLE(III) 

CALL RERAY(MAT1J,AM,FNL,1,IX,2P ) 

IF (KR(17>) 245,265,245 
2^5 CONTTNUF 

WRITF(K0UT,255) 

255 FnRMAT(2X9HDEBUG FLE ) 

WRITF( KCUT,180) (FLE( J) ,J=1 ,MAT1I) 

265 continue 

C determine maximum nonlinear ERRORS 

CALL ABMA X( MAT 1 J, FNLE ,FNLEM, I FNLM) 

ELMM = ABS(ELMM) 

FNLMM=ABS(FNLEM) 

C evaluate nonlinear CORRECTIONS 
DO 615 I=1,maT1J 
615 OVNL( I )=ENL( I ) 

625 CONTINUE 

C evaluate linear CORRECTIONS 
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DO ^30 1 = 1, MATH 
00 f30 J=1,MAT1J 

63C ‘=IF(I) = FL<=(I) - DVNUJ) ♦ BAl(I,J) 
FLF(NF'^^ )=0. 
no 62C T=1,NETA 
L = T-i-META-1 
6?C COP?R( I)=FLE(l) 
continue 

jr f.33 4 CORAPdCORM) / CORMA) 670,675,f75 
67C TF (EA<-F-. 2) 675,671 ,671 

^71 RUMP = 2. ♦BUMP 

675 CALL ABMAXCNFTA , CORAP , CORMA, TCORM) 

IF (KR(17)) 680,680,685 
6FC I'' (KDflO)) 690, “^05, 6^0 
6P5 Continue 

I<R( 17) = KPf 17) - 1 

<*cc comttnuf 

6-=5 FORM/' T( 2X23HDEBUG CPPRECTTCNR FL,NL ) 
WPTTF(KnijT,69'^) 

WPTTE(KCL'T,1®0) (FLE( J) , J=1 ,MAT1 I) 
Wi?I'rF(KrUT,18C) (OVNL( J) ,J = 1,MA'^1J) 

"^0 5 CONTINUE 

C CORRECT FPIMAPY VARIABLES 

nUM=. 1 /RUMO 

EA^F=awTN 1(3.*FA EF ,1. ,0UM/ ABS (CORMfl ) ) 

71C IF (kr(13)) 72C,720,715 
''15 HUM = KP( 13) 

EA?'^ = fiMlNKOLlM / 10.,FASE) 

■'20 IF (=^ASE - 1.0) 725, TAG, ■'AO 
725 00 “'30 1 = 1, math 
73C FLE( I ) = FLE (1 ) * EAEF 
DO 735 I=1,MAT1J 
735 OVNL(T) = DVNL(I) ♦ FA<‘E 
7AG CONTINUE 
■"=0 DO 7a*5 j = 2,NETA 

F( 1,J) = F(l ,J) + PLF< J - 1) 

7 AC rPM 
M = 

on 7P0 j=i,neta 

M = M + 1 

F( 2,J) = F(2,J) + FLF(M - 1) 

78C continue 

on 785 J=2,NETA 
M = M + I 

■'05 F(3,J - 1) = F(3,J - 1) > FlE(m - 1) 

3CC Ff 1,1)=F( l,l)+OVNL(l) 

F( 3,N'^T,« )=F( 3,NF^A)+PVNL(2 ) 

0^ »35 J=1,NETA 
f('*,J) = F(A, J)+DVNL( J^2) 

^35 cnMTTNur 

IF ( I^^S - -^o) °*^C,8A0,850 
P4Q TC (T7-77 - 777) 6A5,P50,RA5 

045 1777 = 777 

TTC = JQ 

85C CONTINU‘D 
RETURN 
END 



SUBROL'-^INE IMONE 
OIMENSTCN RRPI5) 

COMMON /FRPCnv/FLE( 73) ,FLA(73) , F LEM, ELM ( 3 ) , ELMM, IFLM,NE 
1LM,ILMM,0FL( 73) ,FNLE(28) ,ENL(2fi ) , fnLEM , ENLM ( 3 ) , FNLmM, I 
2FNLM,NFNLM,T NLMM,dfnL (28) ,0RNL(6) 

COMMON /FTACOm/eTA (25) ,DE'ta (2A) ,0S0( 2^ ) , OCU ( 2 A ) , B 1 ( 2A ) , 
1B2( 2^) , LAP ( ICO) ,BAI(73 ,28) 

common /HI SCOM/T Z, T 1 ,T2 ,TF(25) ,C 1 ,C2 , C3 , CA , BET a, ZM(5,25 
1) ,XT( 50) ,HF( 25,5) ,H(JF ,HHUE,0LX2 ,(:3M ( 50 ) , B§T a^M ( $0 ) 

COMMON/INTCOM/KR(20) ,kin,kout,matii .matij.neta, I, 15, ns 
I, IT,NTI me ,NAM,MLF 0,NNLF0,NRNL,ITS,CASE( 15) ,B( 3) , MWE,NO 
2N,ITFM,MITFM,KR17 ,NB^.NBT2 , 1 DENT ,UMFS , FREQ, EPS 
CnMMON/NONCOM/AM(2B,28) ,DVNL(2« ) ,C5 , C6 , C7 , C8 , CR , C 10, C 1 
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no 



U,C12 

CnMMDN/PRPC0M/PP,SC.XM(5) 

CnMMON/VARCGM/F(5,25) .FTMl (5 ,25 ) ,FTM2(5,25) 

C EVALUATF groupings which CGNTRIBUTF to n~n PORTION 
C OF COFFFS VARIABLES WITH OIMFNSTON (NFTA-l) 

TF( I-2)4C1,AC02,401 
40C2 XM(5)=C. 

EVALUATE XM, WHICH CONTRIBUTES TO ERRORS AND TO COEPFS 
AT (I) AND ( I-l) 

AOl CALL TAYLOR (OE TA ( I - 1) , F ( 2 , 1 - 1) , F ( 2 , I ) , X M ) 

C EVAL PORTION OF NLE DEPENDENT ON XM,... AND GROUPINGS 

C EVAL at (I-n. 

40 3 CT2 = F( 2,n*XV( l) + F(3 ,n*XM(2)+F(4, I)*XM(3)+F(4t I-1)*XM 

1(4) 

XM( 5)=XM( 5)+C72 

C73=F( 2,n 1 ,T-1)+F(3,I)*ZM(2 n«ZM( 3, I-l) 

1+’=(4,I-1)«ZM(4,I-1 ) 

CNLF( I )=-(-C7-C11+DETA(I-1 )/2.*(0ETA+C^*^5)+.5*C2♦C72- 
12.*C73) 

RRR( 1)=FNLE( I ) 

C EVAL PORTION OF ORIG COEFFS OF AM DEPENDENT UPON PAPAM 

C FVAl AT (I-l). 

AC 5 AM( I , I + 1)=C2*XM(4)-2 .^ZM(/.,i-i) 

4C7 IF( 1-2) 413,417,418 

417 AM( 2,1)=-CP 

C FVAL I-l PORTION OF COEFS, AND FORM CONTRIBUTIONS TO 

C AM AND ERRORS. 

418 LPI=NETA-2+I 
DUM1=-C1C 

421 CALL AMSET (OUMl ,LPI ,2 ,RPR) 

411 LPI=LPI+NETA 
DUM1=-1. 

423 CALL AMSFT ( DUMl ,LPI ,3 ,RPR ) 

IF( 1-2)424, 42B, 424 

424 CALL AMSFT ( -C O , I -2 ,4 , RRR ) 

425 CONTINUE 

RRR( 5)=FNLE( I) 

oqgq FTP MA T ( 2X , I 5 , 2 XI P5 E 1 2 . 4 ) 

IF(»fR(17)) 9501,9902,9901 

9901 CnNTTNUF 
WRTTF(K0UT,9C01) 

ROOl F0 RMAt( 34 H debug FNLE FRQm I MONE PLUS COEFS ) 

WRITE! KOUT, 0999) I ,RRR 

9902 CONTINUE 

DCTypN} 

END 



SUBROUTINE AM$FT ( CA M , LPI , NN , RR R ) 
dimension RRR(5) 

COMMON /EPPC0m/FLE( 73) ,ELA(73) , F LEM, ELM ( 3 ) , ELMM, IFLM,NE 
ILM, ILmm,DFL( 73) ,FNLE (28) ,ENL(2B ) ,fnlFM , ENLM ( 3 ) , ENLMM, I 
2FNLM,NENLM,T NLMM,DFNL(28) ,DRNL(6) 

CnMM6N/FTAC0M/ETA(25) , DETA ( 24 ) , DSO ( 24 ) , DCU ( 24 ) , R 1 ( 24 ) , 
1B2(24) ,LAR(1C0) ,BA1(73,28) 

C0MMnN/INTC0M/KR(20) , KI N ,KOUT , MATl I ,MAT1 J,NEt,A, I , IS, NS 
1, IT,NtimF ,NAM,nlEC,NNLEO,NRNL ,TTS,CASE( 15 ) ,B( 8) , MWE,NO 
2N, ITEm,NI TEM,KR17 ,nBT,NBT2 , 1 DENT ,UMFS , FR FO, »=PS 
C0MM0N/NCNC0M/AM(28,28) ,DVNL (2P ) ,C5 , C6 , C7 , CB , C9 , C 10 , 

1C 11, C 12 

FNLE( ! )=FNLE (I )-CAM*FLF(LPT ) 

IF(NN-l)2,ltl 

1 RRP(NN)=fnlE(I) 

2 no 3 J=1 ,matij 

3 AM( I, J)=AM(I ,j)-cAM«BA1 (LPI ,J) 

R'=T|jPN 

FND 



SUBROUTINE ICNIY 
DIMFNSICN RRR(4) 

C0MM0n/ERRC0M/FLF(73) ,ELA(73) , F lEM , ELM ( 3 ) , ELMM, IFLM,NE 
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c 

c 

c 



c 

r 

r 

c 



lL»^,ILMV,nFL( 73) ,»=NLE ( 29 ) ,ENL { 28 ) ,FNLEM,6MM( 3 ), ENLMM, I 
2»=NLFM,NFNtM, INLMM OFNL (28) |DRNL (6) 

COMMON /5^TAC0V/FTA (25) ,nFTA(24) ,DS0(24) , OCU { 24 ) , B I ( 24 ) , 
18 2(24) ,U*R(100) ,RAir73,29) 

COMM^N/HIEC^^^/'r7,Tl ,T2,TF(25) ,C1,C2 jC 3,C4,8ETA, ZM(5t25 
1 ) , XT ( 93) ,HP( 25 ,5) ,HIF , HHUE , OLX2 , C3 M ( 50 ) .8ETAM (50 ) 
CnMMON/INTCny/KR (20) , KI N ,KCUT , M ATI T , ^AT 1 J . NFT A, I, IS, NS 
1 , ,NA l^,NLFO, NNLEO.NRNLf ITS , CASE ( 1 5 ) , B ( 8 ) , MW = , NO 

2N, tT«=w,MTFM,kR 17,NBT,NBT2 ,I0ENT, ,(IMFS , FPEO, EPS 
CnMM0N/rCNC0M/AM(28t28) , DV NL ( 28 ) ,C5 , C6 , C7 , C8 , CP t C 10 , 



IC11,C12 

COMMCN/PPPCOM/PR ,SC,XM(5) 

rnMMGN/VARCQM/F T*" ,25) ,FTM1 (5 ,25 ) , FT M2 ( 5, 25) 

ADO fPNTR TBUTI ONS OF I '^0 NONLINEAR ERRORS 
evaluate grouping^ which are used only AT I (NOT AT 

FN!^E( I )=FNLE ( I )-(C7 + DETA(I -1) /2 . ’f' ( BET A + C^ *C5 ) +C 1 1 ) 

RRR ( I )=FMLF( I ) 

FVAL OOR'^TON OF GRIG COEFF'^ OF 5m nEPENOENT UPON PARAM 
FVAL at !. 

403 ^M( T,T+2)=AM(I ,I+2)+r2*XM(3)-2.*ZM(3,I-l) 

AOP CONTRIRUTI 0N«^ OF I nEPENOENT CHEFFS TO AM AND 

nonltnfar errors. 

|F( T-NC'ri)43q^^.3p^/^3o 

A?o AM(I , 2 )=AM(I ,2) +1.^C2*XM(2 )-2.*7M(2 ,T-1 ) 

^3P CALL AMSET ( CP , I -1 ,2 , RR R ) 

LDT = n=^ta-1 + I 

0'JM1=C 1C + C2^XM( l)-2.*ZM(l,I-l) 

44A CALL AMSFT ( OUMl , L PI , 3 , RRR ) 

418 IF( I~NFTA )45E,4C‘=,459 
4*5? LR! = ?*NFTA-1 + I 

nUMl=l.+C2*XM( 2 ) «-2.«7v( 2,1 -1 ) 

461 CALL A^SET ( DUMl ,L R! t4 , RRR ) 

405 CPNTINl'P 

990 Q F'^RMAT(2X,I5,2X1P4E12.4) 

IF(KR(17)) PPOl ,P902,9P01 
PPOl CONTINUF 

WRT-F(‘<^U’^,?C01) 

POCl FnpM^T(34H OEPUG ^NLE FROM lONLY PLUS CFEFS ) 

WRJTf (KCUT,99Q9) T,RRR 
P?C2 C'^NTINMIF 



PF'^UPN 

END 



SURRPUTIME OLTdut 

PC*|_ M(t 

OIMFNSTfN C’^(5) ,U( 2*^) ,Y(25 ) 

Cn^MCN/‘^rGCOV/rEM(*^C) ,PTE(5C) ,0UE5 ,UE ( 50 , 50 ) , RE ( 50, ) 

l,RHn,MU,T,nS!P( ■^0) ,IPSIP 

common /FTA r0y/fj^( 25) ,dETA (24) , OSO ( 24 ) , DCU ( 24 ) , B 1 ( 24 ) , 
1B2< 24 ) ,L*R ( 100) ,BA 1( 73 ,2«) 

rnMMQN/HTSC0M/T7,Tl , T2 , TF ( 25 ) , C 1 ,C2 , C3 , C^ , BET A, 7M(5,25 
1) ,XM 50 ,HF(25,5) ,HUF ,HHUE ,0LX2 ,C3 M ( fq ) , R fr am ( 5C ) 

common /INTCQM/kR (2 0) ,KIN,K0UT,MAT1I , mat 1J.NET A, I, IS f NS 
1, IT, NT IMF ,NA M, NLP 0,NNLE0,NRNLi T'^S, CASE ( ) , B ( 8 ) , MW E , NO 

2N, !TFM,mtfM,KP17 ,NBt,NBT2 , I OENT ,UMFS , FR fq, fos 

CnMM0N/N0NCnM/AM(28,28) ,DVNL(28 ) ,C5 , C6 , C? , C8 , CP , C 10 , 
1C11,C12 

CnMMGN/FRMCOM/Ti me (^C) ,PRF(50), RTET , S ( ) , RCK A^^ ( 50 ) , RN 
lOSF , VKAR ,NOT SC , lOT SC (50) ,NSO(10 ) ,MSD( 10 ) , lopF^PADNO, CO 
2NP,DEL'^, CYCLE 

C0MM0N/PRPCCM/PR,SC ,XM(E ) 

CnMMGN/TFMC0M/SDDUM2 (10) , DPR (50) , OUMMI ( 25 ) , S LOPE ( 25 ) , R 
IFOUMf 2^) ,S0UV1 ( =^0) ,SnUM2 ( 50) ,FW0UM(50 ) ,X ICON (50 ), FWCO 
2N( *"0) ,nunS(5C) 

COMM0N/VARC0m/F( 5 ,25) ,FTM1 ( 5 ,25 ) ,FtM2 ( 5 , 25 ) 

C0MM0N/WALC0M/FW(50) ,RHOVW(50) , TFW, IRHOVW, JRHGVW 
1 F0PMAT( /<»X65HFDGE STATIC XI 

1 flux norm ROKAP /6X57HVEL0CITY PRESSURE 

2(LR/ PPTA parameter /6X47H( FT/S*=C ) 

3(ATM) SEC)*«2 (FT) 
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WALL 



3C0 

3C2 

3C1 



3011 

3012 
3051 



1*501 

190 

303 

3CA 

182 



212 

214 



183 

1930 

194 



4 /3X6( 2X1PE1C.3) /) 

FORMAT( /10X54HWALL FR. MQH TRANS 

1 mass FLUXFS /11X56HCOEFF, COEFF, SHEAR 

2 meCH total /13X53HCF RHO*UE«CF/2 

3 REM GAS /35X31H(IB/SQ FT) 

4 (LP/SO FT SEC) /6X5 ( 3 XIPEIO .3 ) / ) 

forma T( /10X58HREYN0LDS STROUHAL MOMENTUM D 

IISOLACF. SHAPE /I 1 X58HNUMBER NUMBER THIC 

2KNFSS, THICKNESS, FACTOR, /10X60H 

3 theta OFLSTAP DELSTAR/ /38X30H(F 

4T) (FT) THETA /4X5(3X19E10.3)/ ) 

F0RMAT( /3l X17HNODAL INFORMATION // 3X50HOISTANCE 

1 STREAM VELOCITY SHEAR /2X75HFRQM WALL 

2 ETA FUNCTION RATIO FUNCTION ^PPP 

3 VELOCITY /5X72H(FT) (F) (FP=U/UE 

4) IFF9) (FT/SEC) /) 

FORMAT! 1X7(IX1PE10.3) ) 
cOPMAT( 1FPFIC,3) 

DATA IPLANK/2H / 

IF (KR(ll)-l) 300,300,301 
(KR(C))3Q2j3C2,301 
OHOVWI IS)=C1^F(1,1)^'HF(1 ,5) 

C R9=-C 3M( I *;) «MU 
ouMi=-i./c3Mn$) 

01JM2 = RH0VW( I S) /C3M( I S) 

IF ( UE( IS, ITEM)- 1. ) 3012, 3011 ,3012 
UE( TS,TTFM)=r. 

CnNTTNUF 

continue 

OEP( 1)=C.C 
OEP( 2)=C.C 
VMECH=-PUM2 
OUw^=VMFCH*lCO, 

IF (OUM4-DUM2) 1901,190,190 
VMFCH=0. 

TF( ABS(BETA)-.OOCI) 303,304,30^' 

BFTA=C. 

Y( 1)=0. 

DO 182 T=2,NETA 

Y( I )=Y( I-1)+C89«0ETA(I-1) 

SHFAP=MU«UEI IS, ITEM) /C89 ♦F ( 3 , 1 ) /32 . 17405 

WOTTF(KGUT,l)UEnS,ITEM) ,PF(fs,ITEM) , X I (I S ) , B ET A , OUM 1 , 
IROKAO( TS) 

nUMl=RHO«UE( IS,I^EM) /My 

CC=Mn/CP9+F(3,l) 

DUM4^1. 

RFS=nUMl«S( I S) 

WCC=CF«2. /RHG/UE (I 5,!TEM) 

0UM3=Co 9«(F( l,NF-rA)-F(l ,1) ) 

0ELST=V(NETA)-DUM3 

redfls=dumi«dflst 

THM0M=dUV3-C 69*XM( 5) 

R cTHMn=nijMl^THMOM 

shapf=oelst/thwcm 
STR=freo/UE( IS, ITEM) 

DUMl = OUMi=«tS( TS) 

STR=ftr«S(TS) 

WR!TF( KCUT,2 ) WFC ,CF , SHEAR , VMECH , DUM2 
WP T TF( kput,3)DUM1 ,STR ,THMQm, oeL?^T , shape 
WPITF(kCUT,4) 

DO 183 I=1,NFTA 

COMPUTE true VALUES OF F(I,J) AND ETA 
REF1=P( 1,1 ) 

9EF2=F( 2,1 ) 

PEF3=F( 3,1 ) 
rcf4=F(4,T) 

PFTA=FTA( I ) 

U(!)=PEF2«UE(IS,ITFM) 
ctan=RETA /(FTA (NETA) ) 

WRITE ( KCUT,5)Y(I ) , RET A , REF 1 , REF 2 , REF 3 , R EF4 , U ( I) 

IF ( IDENT-TRLANK) 194,309,194 
I9(N0N) 3C9,308,3C9 
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3CP nUM3 = RH0*UEnS,TTEW) 

DU'^A=-1, /C3 
3C9 CONTIMIT 

I^f KP( 1^) )31C,311.310 

3 1C WRTTF(KrUT,^)( (F(! ,.J) tJ = l»NFTA) tl = l,4) 
311 R'=TIJPK! 

RNO 
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